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ABSTRACT 


This  thesis  formulates  the  full  nonlinear  equations  of  motion  for  determining  the 
stability  of  helicopter  coupled  rotor-fuselage  motion  utilizing  MATLAB®‘s  Symbolic 
Math  Toolbox.  Using  the  extended  symbolic  processor  toolbox,  the  goal  of  this  work  was 
to  eliminate  the  time  consuming  process  of  converting  Fortran  or  C  code  generated  by  the 
symbolic  processor,  MAPLE®  into  a  MATLAB®  useable  format  where  it  is  further 
incorporated  into  an  ‘S-function’  to  be  used  in  the  dynamic  simulation  environment. 

The  formulation  of  the  equations  of  motion  utilized  in  this  process  is  unique  in 
that  it  uses  the  complete  set  of  nonlinear  terms  in  the  equations  of  motions  without 
utilizing  ordering  schemes,  small  angle  assumptions,  linearizing  techniques,  or  other 
simplifying  assumptions.  After  derivation,  the  equations  of  motion  are  numerically 
integrated  using  the  dynamic  simulation  software  SIMULINK®  and  a  time  history  plot  is 
generated  of  blade  and  fuselage  motion.  The  equations  of  motion  are  regenerated  with 
each  time  step  allowing  the  adjustment  of  characteristic  structural,  blade  and  dampening 
properties.  These  time  traces  can  be  used  to  explore  the  effects  of  damping 
nonlinearities,  structural  nonlinearities,  active  control,  individual  blade  control,  and 
damper  failure  on  ground  resonance. 
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I. 


INTRODUCTION 


A.  DISCUSSION 

Air  and  ground  resonance  are  potentially  destructive  mechanical  instabilities  that 
may  occur  in  helicopters  due  to  dynamically  coupled  interaction  between  the  rigid  body 
motion  of  the  fuselage  and  oscillations  of  the  rotor  blades  in  their  plane  of  rotation.  For 
over  40  years  engineers  have  felt  reasonably  comfortable  with  their  ability  to  predict  and 
control  air  and  ground  mechanical  instabilities.  The  self  excited  vibrations  that  occur 
while  an  aircraft  is  in  contact  with  the  ground  is  known  as  ground  resonance  and  this  will 
be  the  primary  focus  of  this  thesis. 

Helicopter  manufacturers  employ  several  methods  of  blade  root  retention  to 
obtain  the  flapping  and  lead-lag  degrees  of  freedom  necessary  for  a  helicopter  rotor  blade 
to  operate  properly.  These  methods  include:  (1)  fully  articulated;  (2)  bearingless;  (3)  and 
hingeless  main  rotor  designs. 

Fully  articulated  rotor  heads  utilize  two  hinges  that  rotate  with  and  are  attached  to 
each  rotor  blade  in  order  to  provide  the  necessary  freedom  of  movement  for  the  dynamic 
motion  of  the  rotor  blade.  One  hinge  is  positioned  in  a  horizontal  plane  and  allows  for 
flapping  motion  of  the  rotor  blade.  The  other  hinge  is  positioned  in  a  vertical  plane.  This 
hinge  allows  for  lead-lag,  movement  of  the  blade  in  its  plane  of  rotation.  The  lead-lag 
hinge  generally  incorporates  some  damping.  This  is  generally  a  hydraulic  air-oil  strut  or 
an  elastomeric  design,  to  constrain  the  movement  of  the  individual  rotor  blades.  Early 
designs  of  these  dampers  have  tended  to  be  heavy  and  maintenance  intensive. 
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The  hingeless/bearingless  rotor  is  a  variation  of  the  fully  articulated  rotor  system. 
In  this  configuration  the  designer  has  substituted  a  flexible  section  to  replace  both  the 
lead-lag  hinge  as  well  as  the  flapping  hinge. 

In  short,  by  eliminating  the  need  for  hinges,  improvements  in  technology, 
particularly  composite  materials,  have  allowed  designers  to  significantly  reduce  the 
weight  and  complexity  of  modem  helicopter  rotor  systems.  However,  this 
reduction/elimination  of  hinges  has  also  reduced  the  likely  attachment  points  for  linear 
dampers  and  other  means  of  controlling  the  motion  of  the  rotor  blades.  This  further 
complicates  eliminating  ground  resonance,  as  well  as  the  associated  complexity  of 
manufacturing  and  maintenance  of  the  dampers.  By  simplifying  and  facilitating  the 
process  of  modeling  ground  resonance  the  hope  is  to  open  the  study  of  ground  resonance 
to  a  deeper  understanding  as  well  as  investigate  methods  of  controlling  it.  By  allowing 
the  variation  of  parameters  and  by  receiving  immediate  feedback  as  to  the  resultant 
effects,  this  program  and  future  versions  will  help  improve  the  properties  of  existing 
damper  systems,  facilitate  the  possibility  for  application  of  nonlinear  dampers,  potentially 
eliminate  external  dampers  using  existing  material  technology,  and  facilitate  the 
investigation  of  other  aero-mechanical  instabilities. 

B.  BACKGROUND 

In  the  early  1940’s  the  phenomenon  of  ground  resonance  was  already  being 
investigated  by  NACA.  Robert  Coleman’s  pioneering  work  was  the  first  definitive 
analysis  to  correctly  address  the  issue  of  ground  resonance  [Ref.  I]  and  [Ref.  2].  In  his 
work  Coleman  identified  ground  resonance  to  be  a  self-excited,  mechanical  instability 
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phenomenon.  He  found  the  primary  modes  to  be  excited  in  the  normal  operation  of  the 
rotorcraft  were  the  hinged  rigid  body  response  of  the  blades,  in  their  plane  of  rotation 
coupled  with  horizontal  deflection  of  the  main  rotor  pylon. 

M.  L.  Deutsch  simplified  the  results  of  Coleman  and  Feingold  for  the  practicing 
engineer.  Based  on  Coleman’s  theory  he  developed  a  quantitative  analysis  of  the 
damping  required  to  keep  the  helicopter  free  of  the  instability,  often  referred  to  as 
Deutsch’ s  Criteria  [Ref.  3]. 

Coleman  and  Feingold’ s  work  became  the  basis  for  the  evolution  of  theory  and 
design  techniques  used  for  dealing  with  ground  resonance.  Although  the  classic  theory 
offers  much  insight  and  understanding  into  the  phenomenon,  especially  for  conventional 
articulated  rotor  systems,  the  increasing  popularity  of  hingeless  and  bearingless  rotor 
designs  in  modem  helicopters  called  for  increasingly  more  sophisticated  analytical  tools. 

More  sophisticated  analytical  tools  became  possible  with  the  evolution  of  digital 
computers  and  the  improvement  in  computational  power.  Peters  and  Hohenemser  [Ref. 
24]  apply  Floquet  analysis  to  the  problem  of  lifting  rotor  stability.  Floquet  analysis  is  a 
method  which  can  be  used  to  determine  the  stability  of  solutions  to  systems  of  linear 
ordinary  differential  equations  with  periodic  coefficients.  The  Floquet  transition  matrix 
which  relates  the  system  state  variables  at  the  beginning  and  end  of  a  rotational  period  is 
computed  by  numerical  time  wise  integration.  The  eigenvalues  of  the  transition  matrix 
are  a  measure  of  system  stability.  Hammond  [Ref.  25]  applies  Floquet  analysis  to  the 
prediction  of  mechanical  instabilities,  specifically  examining  the  case  of  unbalanced  lead- 
lag  damping.  The  unbalanced  problem  requires  solution  of  the  equations  of  motion  with 
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the  blade  equations  expressed  in  the  rotating  reference  frame  because  a  transformation  to 
the  fixed  system  is  no  longer  possible  for  a  ground  resonance  analysis  as  was  possible  for 
the  isometric  case.  As  a  result,  you  are  left  with  a  system  of  equations  with  periodic 
coefficients  which  can  be  handled  by  the  Floquet  method. 

Hingeless  and  bearingless  rotor  configurations  often  face  the  additional  difficulty 
of  air  resonance.  Aerodynamics  may  play  more  than  a  passive  roll  in  the  ground 
resonance  regime  in  hingeless  systems  in  contrast  to  articulated  systems  where 
aerodynamics  have  little  effect.  As  a  result,  more  complex  models  are  required  to 
accurately  represent  the  physics  of  the  helicopter  aeromechanical  stability  problem. 
Models  must  include  blade  flap  and  torsional  degrees  of  freedom  as  well  as  lead-lag 
degrees  of  freedom.  Fuselage  models  also  should  include  pitch  and  roll  as  well  as 
translational  degrees  of  freedom.  Aerodynamic  models  can  range  from  quasi-steady  strip 
theory  to  unsteady  aerodynamic  theories  which  include  elaborate  wake  models  or 
dynamic  inflow  models.  Ormiston  [Ref.  26]  utilizes  a  rigid  blade  and  rigid  fuselage 
model  with  flap-lag  and  pitch-roll  degrees  of  freedom  to  conduct  parametric 
investigations  based  on  an  eigenvalue  analysis.  As  is  typical,  the  equations  of  motion 
were  derived  by  a  Newtonian  approach  and  the  resulting  system  of  nonlinear  differential 
equations  are  linearized  for  small  perturbations.  The  model  includes  linear  rotor  blade 
and  landing  gear  springs,  viscous  damping,  and  quasi-steady  aerodynamics.  Freidmann 
and  Venkatesan  [Ref.  27]  and  Freidmann  and  Warmbrodt  [Ref.  14]  derive  the  complete 
set  of  governing  equations  of  a  helicopter  rotor  coupled  to  a  rigid  body  fuselage.  The 
equations  account  for  rotor  blade  elastic  deformations  and  include  quasi-steady 
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aerodynamics  or  modified  Theordorsen  unsteady  aerodynamic  theory.  In  deriving  the 
full  equations  of  motion,  Freidmann  et  al.,  stress  the  importance  of  applying  an  ordering 
scheme  to  the  process  in  order  to  handle  the  complexity  of  the  equations  and  enormous 
number  of  terms  generated  by  their  expansion.  The  equations,  as  presented  by  Freidmann 
et  al.  [Ref.  27  and  14],  are  in  a  form  which  makes  them  generally  applicable  to  a  wide 
range  of  rotorcraft  problems. 

Another  interest  in  the  study  of  helicopter  ground  resonance  is  the  effect  that 
nonlinear  elastic  and  damping  forces  have  on  stability.  Tongue  [Ref.  23],  Tongue  and 
Flowers  [Ref.  9  and  8],  Tongue  and  Jankowski  [Ref.  28],  and  Tang  and  Dowell  [Ref.  29], 
use  variations  of  the  nonlinear  technique  of  harmonic  balance  using  describing  functions 
to  represent  nonlinear  damping.  The  technique  is  useful  for  investigating  limit  cycle 
behavior  of  strongly  nonlinear  systems  and  its  impact  on  system  stability. 

Helicopter  aeromechanical  instabilities  can  be  analyzed  by  methods  ranging  from 
Coleman’s  classic  analysis  to  direct  time  integration  of  the  equations  of  motion.  As  . 
engineers  strive  to  develop  rotor  systems  free  of  ground  and  air  resonance  which  do  not 
require  the  addition  of  maintenance  intensive  mechanical  damping  systems,  more 
elaborate  models  will  be  needed  to  accurately  capture  all  physical  aspects  of  the  problem. 
To  achieve  the  truly  damperless  rotor  Ormiston  [Ref.  30]  addresses  three  different 
approaches  which  may  be  feasible,  1)  incorporating  high  damping  material  into  the  blade 
or  flexbeam  structure,  2)  automatic  feedback  control,  and  3)  development  of  aeroelastic 
couplings  to  provide  inherent  stability.  These  three  approaches  provided  the  impetus 
behind  the  work  performed  by  LT  Christopher  S.  Robinson. 
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A  full  nonlinear  simulation  model  of  coupled  rotor/fuselage  interaction  was 
created  by  Robinson  in  March  of  1997  at  the  Naval  Postgraduate  School  (NPS)  [Ref.  4]. 
This  model  is  a  dedicated  Coleman  analysis  tool  that  initially  utilizes  a  symbolic 
processor,  Maple  to  derive  the  full  nonlinear  equations  of  motion  of  a  helicopter  using 
the  LaGrange  equation.  In  Robinson’s  thesis  the  derived  equations  of  motion  are  then 
converted,  in  Maple®,  to  Fortran  or  C  and  then  converted  into  a  MATLAB® 
programming  language.  The  converted  MATLAB®  result  is  incorporated  into  a 
SIMULINK  S-function,  producing  time  history  plots  of  blade/fuselage  motion.  This 
process  is  unique  in  that  it  uses  the  complete  set  of  nonlinear  terms  in  the  equations  of 
motions  without  utilizing  ordering  schemes,  small  angle  assumptions,  linearizing 
techniques,  or  other  simplifying  assumptions. 

Further  development  of  the  NPS  modeler  was  accomplished  by  Robinson  with  the 
help  of  LCDR  Robert  L.  King,  [Ref.  12, 18, 19,  20].  The  modeler  was  used  to  simulate 
the  Froude  Scale  RAH-66  [Ref.  21].  It  was  also  used  to  investigate  an  interblade 
coupling  approach  to  improved  rotor  stability  without  lag  dampers  [Ref.  22]. 

King  completed  his  dissertation  in  June  1999  [Ref.  5].  In  his  dissertation.  King 
adopted  Robinson’s  work  to  further  explore  the  potential  of  eliminating  the  snubber- 
damper  or  damper  on  hingeless  rotor  designs  and  replace  it  with  a  flexbeam  that  has  been 
modified  to  possess  nonlinear  properties.  It  is  King’s  research  in  this  field  that  has 
provided  the  motivation  to  accomplish  the  following  work. 

Since  Robinson’s  original  work,  MATLAB®  has  incorporated  a  symbolic 
processing  toolbox  into  its  software.  By  using  the  extended  symbolic  processor  toolbox 
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in  MATLAB®  the  goal  is  to  eliminate  the  awkward  and  time  consuming  process  of 
converting  the  Fortran  or  C  code  generated  by  MAPLE®  into  a  MATLAB®  useable 
format  where  it  is  further  incorporated  into  an  ‘S-function’  to  be  used  in  the  dynamic 
simulation  environment.  Direct  interface  between  the  derivation  of  the  equations  of 
motion  and  the  actual  simulation  process  results  in  fewer  steps  required  to  develop  a 
system,  it  reduces  the  computer  programs  required  to  run  a  simulation,  and  reduces  the 
potential  contamination  of  the  solution  that  may  result  from  human  interaction. 

Hopefully  this  will  encourage  further  use  of  the  rotor-fuselage  coupled  motion  simulation 
and  enable  researchers  of  all  levels  to  explore  the  phenomenon  of  ground  resonance. 

After  a  brief  explanation  of  the  work  performed  by  Robinson,  Rafanello,  and 
King  an  introduction  to  MATLAB®,  MATLAB®  the  symbolic  processing  toolbox,  the  S- 
function,  and  the  Simulink  model  will  be  provided.  The  goal,  again,  is  to  provide  a  path 
for  future  use  of  the  NPS  modeler. 
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II.  CURRENT  MODEL 


A.  REVIEW  OF  WORK  PERFORMED  BY  ROBINSON 

1.  Simple  Model 

As  with  this  thesis,  Robinson’s  work  began  with  a  simplified  model  using  reduced 
degrees  of  freedom  on  a  three  bladed  model.  The  simplified  model  is  based  on  that  of 
Coleman  as  is  shown  in  Figure  1,  [Ref.  1].  Elastic  forces  generated  by  rotor  blade  and 
flexbeam  motion  were  modeled  as  a  linear  torsional  spring  located  at  the  effective  hinge 
position  of  the  blade.  Landing  gear  stiffness  was  also  modeled  with  linear  springs.  The 
landing  gear  and  lead-lag  dampers  were  modeled  with  linear  dashpot  type  dampers. 
Transformation  between  the  various  systems,  in  order  to  develop  Lagrange’s  equations  of 
motion,  were  accomplished  using  Euler  angle  rotations.  For  the  simple  model  the 
coordinate  transformations  used  were  hub  (parallel  to  fuselage  but  offset  a  distance  h  in 
the  positive  z  direction)  to  inertial  (fixed  relative  to  the  earth),  blade  undeformed  (fixed  to 
the  effective  hinge  position)  to  hub,  and  blade  deformed  (fixed  to  the  effective  hinge 
position  with  the  x-axis  coincident  with  the  blade  ‘elastic’  axis)  to  blade  undeformed. 
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Figure  1.  Simplified  Rotor  Model  [From  Ref.  4] 


Robinson  showed  excellent  agreement  between  his  model  and  the  Coleman 
model,  with  departure  occurring  only  when  displacements  are  very  large.  Figure  2.  This 
is  to  be  expected  since  Robinson’s  model  does  not  assume  small  angle  theory  whereas  the 
Coleman-Feingold-Bramwell  model  does. 
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Figure  2.  Validation  of  model  with  solution  of  Coleman’s  model  [From  Ref.  4] 


Figure  3  is  a  parametric  plot  of  the  simple  model  with  an  isotropic  pylon  and 
rotor.  The  damping  ratio  from  the  moving  block  result  is  plotted  versus  oo/cof  for  various 
Deutsch  numbers  [Ref.  15]. 
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Figure  3.  Validation  of  model  against  Deutsch  Criteria  [After  Ref.  4] 

This  model  allows  for  the  following  degrees  of  freedom: 
ul  -  Fuselage  translation  in  1-direction  (x-direction) 
u2  -  Fuselage  translation  in  the  2-direction  (y-direction) 

Ck  -  Lead-lag  angular  displacement  of  k*  rotor  blade 

2.  Complex  Model 

The  complex  model  is  based  on  that  used  by  Straub  [Ref.  6].  This  model  assumes 
rigid  blades  and  fuselage,  as  in  the  simple  model,  with  flap  and  lead-lag  hinges  as 
coincident.  It  incorporates  fuselage  rotation  as  well  as  blade  flap,  requiring  additional 
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transformations.  Figures  4  and  5  show  a  schematic  representation  of  the  complex  model 
and  the  transformations  utilized  in  the  derivation.  The  transformations  required  in  the 
complex  model  are  fuselage  (fixed  to  center  of  gravity  of  fuselage)  to  inertial,  hub  to 
fuselage,  blade  undeformed  to  hub,  blade  deformed  to  blade  undeformed. 
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Figure  5.  Schematic  representation  of  complex  model  transformations 
This  model  utilizes  the  following  degrees  of  freedom: 

ul  -  Fuselage  translation  in  1 -direction  (x-direction) 
u2  -  Fuselage  translation  in  the  2-direction  (y-direction) 
rl  -  Fuselage  rotation  about  1-axis  (roll) 
r2  -  Fuselage  rotation  about  2-axis  (pitch) 

£k  -  Lead-lag  angular  displacement  of  kth  rotor  blade 
Pk  -  Flap  angular  displacement  of  k*11  rotor  blade 
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3.  Derivation  of  the  Lagrange  Equation 

A  brief  overview  of  the  derivation  accomplished  by  Robinson  will  be  included  in 
this  thesis.  For  a  more  extensive  explanation  see  Robinson’s  thesis  [Ref.  4].  The 
Lagrangian  equation  may  be  expressed  as  follows: 

d  f  dT  "j  dT  |  dU  |  dD  _ 
dt{dq,  )  d<?  +  dq,  +  dqt 

Where  T  is  the  kinetic  energy,  U  is  the  potential  energy,  D  is  the  dissipation 
function,  Fj  is  the  generalized  force,  and  q*  is  the  generalized  displacement.  In  the 
complex  model  the  generalized  force  term  Fj,  will  describe  the  aerodynamic  forces  on  the 
individual  rotor  blades.  This  derivation  develops  a  system  of  homogeneous  equations. 
The  various  energy  terms  are  broken  down  into  two  categories,  blade  motion  and 
fuselage  motion  to  give  the  following  equations: 

r=rr+£(r„). 

*=1 

v  =u  F +£<{/  A 
*=1 

D  =  Df+  £(£/,), 

k=  1 

Where  the  subscripts  F  and  B  indicate  the  fuselage  and  rotor  blade  respectively. 
a.  Kinetic  Energy  Terms 

The  kinetic  energy  of  the  kth  rotor  blade  is  given  by  the  following 

expression: 
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R  1 

C TB)  =  \-m\p-p)dr 
0  2 

Here  p  is  the  position  of  a  point  on  the  elastic  axis  of  the  kth  rotor  blade 
with  respect  to  the  inertial  coordinate  system  at  any  instant  in  time,  and  m’  is  the  mass 
distribution  per  unit  length  of  the  blade  (mass  distribution  per  unit  length  is  assumed  to 
be  uniform).  The  position  of  a  point  on  the  elastic  axis  of  a  rotor  blade,  p ,  is  expressed 
as  the  sum  of  relative  positions  with  respect  to  the  various  coordinate  systems 
transformed  to  the  inertial  system.  Thus, 

P  =  (Pr_'),  +(Ph_F^1  +(Pbu_h  \  +(.PBd_Bu\  +(Pp_BD^n 

where,  for  example,  the  term  (pBu  H ),  is  the  position  of  the  origin  of  the 

undeformed  blade  coordinate  System  with  respect  to  the  hub  coordinate  system 
transformed  into  the  inertial  coordinate  system. 

Following  the  proper  transformations,  the  resultant  expression  provides 
the  position  of  an  arbitrary  point  on  the  elastic  axis  of  the  kth  rotor  blade,  with  respect  to 
the  inertial  coordinate  system,  at  any  instant  in  time,  in  terms  of  the  system  degrees  of 

freedom.  The  time  derivative  of  this  expression  gives  the  velocity,  p ,  which  is 
substituted  into  the  equation  for  kinetic  energy  for  the  kth  rotor  blade. 

The  fuselage  kinetic  energy  in  terms  of  translational  and  rotational  degrees 

of  freedom  is 

OV)™, 

) rot  =-Vl  ^  ^22 ^2  — 

b.  Potential  Energy  Terms 

The  elastic  forces  generated  by  rotor  blade  motion  give  rise  to  a  potential 
energy  term  in  the  Lagrange  equation.  Since  a  rigid  blade  model  was  assumed,  the 
potential  energy  was  modeled  using  equivalent  torsional  springs  to  restrain  the  rotor 
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blade,  with  spring  constants  selected  to  approximate  elastic  forces  due  to  in  plane  and  out 
of  plane  bending  of  the  rotor  blade  (and  the  flexbeam  in  the  hingeless  case).  The 
potential  energy  of  the  kth  rotor  blade  is 

The  fuselage  potential  energy  in  terms  of  translational  and  rotation  degrees 

of  freedom  is 

An  explanation  of  the  validity  of  using  an  equivalent  torsional  spring 
system  to  model  the  elastic  forces  of  a  deformed  rotor  blade  is  given  in  some  detail  in 
Venkatesan  and  Friedmann  [Ref.  7], 

c.  Dissipation  Function  Terms 

System  damping  is  modeled  in  energy  form  by  use  of  a  dissipation 
function,  which  for  the  kth  rotor  blade  of  he  complex  rotor  model  is 

The  dissipation  function  for  the  fuselage  in  terms  of  translational  and  rotational 
degrees  of  freedom  is 

(DF)„„=^Ctuf+±C2il 
(Dr)„=|cR,^+icR2  r? 

d.  Resultant  Equations 

It  is  important  to  note  that  when  applying  Lagrange’s  equation  in 
MAPLE®,  derivatives  with  respect  to  the  degrees  of  freedom  and  the  time  rates  of  change 
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of  the  degrees  of  freedom,  the  time  functional  notation  which  represents  these  variables 
must  be  converted  to  independent  variable  notation.  For  example,  the  flap  angle  degree 


of  freedom,  |3(t),  and  its  time  rate  of  change. 


mo 

dt 


,  would  have  to  be  replaced  in  all  of 


the  energy  expressions  by  the  independent  variables,  |3  and  d(3  respectively,  in  order  for 


ar 


ar 


,  Ml 

Hi 


terms  such  as  and  ^  (where  q;  =  P(t)  and  dt  )  to  be  evaluated  properly  by 


rdT ^ 


,to 


the  MAPLE®  symbolic  engine.  Additionally,  for  the  time  derivative  term,  ^ ' 
be  evaluated  properly,  all  degrees  of  freedom  expressed  in  independent  notation  must  be 
converted  back  to  time  dependent  notation. 

The  equations  of  motion  generated  for  the  simplified  model  (Coleman 
model)  in  this  format  were  compared,  by  Robinson,  to  the  equations  used  by  Flowers  and 
Tongue  [Ref.  8,9].  The  equations  were  found  to  match  exactly  except  for  the  nonlinear 
damping  terms  that  were  not  included  initially  but  are  later  included  in  the  simulation. 
Robinson’s  work  has  been  further  utilized  by  LT  Salvator  Rafanello  and  LCDR  Robert  L. 


King. 


B.  REVIEW  OF  WORK  PERFORMED  BY  RAFANELLO 

Rafanello’ s  work  further  validated  the  NPS  ground  resonance  modeler  by 
matching  simulation  results  with  Professor  Roberts  E.  Wood’s  HSS-2  ground  resonance 
analysis  [Ref.  10, 11].  Rafanello  used  a  five  bladed  version  of  the  simplified  model 
adapted  to  the  H-3  Sea  King.  A  mass-spring-damper  system  of  the  H-3’s  landing  gear 
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was  examined  to  obtain  parameters  to  enter  into  the  NPS  simulation.  A  thorough 
analysis  of  the  natural  frequencies  in  the  lateral  and  roll  modes  at  three  power  settings  of 
0/20/80  percent  airborne  were  made  with  respect  to  the  literature  prepared  by  Coleman, 
Feingold,  and  Deutsch.  Stability  charts  were  developed  from  the  tailored  data  of  the  Sea 
King  and  compared  with  the  modeler’s  results.  Through  these  results  Rafanello  was  able 
to  directly  link  the  classical  work  with  the  NPS  modeler  in  the  roll  and  lateral  modes  of 
the  H-3. 

In  his  work,  Rafanello  was  also  able  to  verify  the  conservative  nature  of 
Deutsch’s  Criteria  as  was  presented  in  Robinson’s  thesis.  He  was  able  to  determine  that 
Deutsch’s  Criteria  is  conservative  in  that  the  critical  point  for  his  criteria,  or  bucket  of 
each  curve,  still  shows  positive  damping  when  analyzed  with  the  NPS  simulation.  Please 
refer  to  Reference  10  for  further  analysis. 

C.  REVIEW  OF  WORK  PERFORMED  BY  KING 

King’s  work  focused  on  examining  alternative  designs  for  a  damperless  helicopter 
rotor.  He  explored  the  potential  of  eliminating  the  snubber-damper  or  damper  on 
hingeless  rotor  designs  and  replacing  it  with  a  flexbeam  that  is  modified  to  possess 
nonlinear  properties.  The  NPS  ground  resonance  modeler  is  well  suited  to  this  task  in 
that  it  can  accurately  model  nonlinear  mechanical  properties.  King  included  nonlinear 
properties  in  the  blade  root  to  produce  potentially  acceptable  bounded  responses  in  the 
parameter  region,  with  the  NPS  modeler,  where  linear  theory  would  have  predicted 
instability.  He  was  also  able  to  determine  regions  of  unacceptable  response  where  linear 
theory  would  have  predicted  stability. 
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King  also  looked  at  linearly  linked  rotor  systems  and  unevenly  spaced  rotors 
using  the  NPS  modeler.  By  linearly  linking  the  rotors  he  was  able  to  change  the  mode 
shapes  of  the  rotor  blades’  lead/lag  motion.  This  alters  the  blade  natural  frequency  in  the 
lead/lag  and  regressing  mode  of  the  rotor.  Therefore,  by  linking  the  blades,  the  potential 
exists  to  avoid  ground/air  resonance  by  detuning  the  rotor-body  dynamics.  This  method 
avoids  ground  and  air  resonance  in  a  similar  manner  to  the  nonlinear  stiffness  approach. 
Unevenly  spacing  the  rotor  blades  produced  the  same  mechanical  instabilities  as  were 
produced  without  altering  the  rotor  configuration,  for  all  cases  tested. 

King  was  further  able  to  use  the  NPS  modeler  to  simulate  the  performance  of  a 
1/6  scale  Comanche  rotor  system.  The  goal  of  the  simulation  was  to  match  the  fixed 
system  damping  values  to  the  test  data  and  to  UMARC  (computer  simulation  modeler 
designed  at  University  of  Maryland).  In  this  comparison  the  NPS  modeler  shows  slightly 
better  results  than  UMARC  at  the  lower  RPM’s  tested.  At  higher  rotor  speed,  the 
difference  between  the  two  simulations  is  negligible.  Figure  6. 
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Figure  6.  Final  NPS  simulation  results  with  increased  geometric  gain  [From  Ref.  5] 


From  these  results  King  concluded  that  “the  NPS  modeler  shows  potential  in 


modeling  rotors  with  nonlinear  physical  parameters  and  further  development  should  be 
considered”  [Ref.  5]. 
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III.  MATLAB® 

A.  INTRODUCTION  AND  HELP 

The  NPS  modeler  has  been  shown  to  be  effective  in  modeling  coupled  rotor- 
fuselage  motion.  In  order  to  help  facilitate  the  functionality  of  the  simulation  program, 
the  inclusion  of  a  symbolic  processing  toolbox  in  MATLAB  provides  a  means  of  refining 
the  simulation  process.  The  computer  programming  contained  within  this  thesis  is 
accomplished  entirely  in  the  computer  programs  MATLAB®  and  SIMULINK®.  Both 
programs  are  developed  by  Mathworks,  Inc.  and  are  run  interchangeably.  MATLAB®  is 
an  array  and  matrix  based  program  allowing  programming  features  similar  to  other 
computer  programming  languages.  In  addition  to  its  matrix  orientation,  MATLAB®  also 
offers  Graphical  User  Interface  (GUI)  tools  and  an  excellent  graphics  package. 
MATLAB®  stores  all  data  as  arrays,  which  lends  itself  to  the  manipulation  of  sets  of  data 
in  a  wide  variety  of  ways.  Recently  Mathworks,  Inc.  has  incorporated  the  Symbolic  Math 
Toolboxes  into  MATLAB®’ s  numeric  environment.  These  toolboxes  supplement 
MATLAB®’ s  numeric  and  graphical  facilities  with  several  other  types  of  mathematical 
computation: 


Facility 

Covers 

Calculus 

Differentiation,  integration,  limits,  summation,  and  Taylor 
series 

Linear  Algebra 

Inverses,  determinants,  eigenvalues,  singular  value 
decomposition,  and  canonical  forms  of  symbolic  matrices 

Simplification 

Methods  of  simplifying  algebraic  expressions 

E'V"; 

Symbolic  and  numerical  solutions  to  algebraic  and 

differential  equations 

Solution  of 

Symbolic  and  numerical  solutions  to  algebraic  and 

Equations 

differential  equations 

Variable-Precision 

Numerical  evaluation  of  mathematical  expressions  to  any 
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Arithmetic 

specified  accuracy 

Transforms 

Fourier,  Laplace,  z-transform,  and  corresponding  inverse 
transforms 

Special  Mathematical 
Functions 

Special  functions  of  classical  applied  mathematics 

Table  1.  MATLAB R  Symbolic  Math  Toolboxes  [After  Ref.  13] 


The  computational  engine  underlying  the  toolboxes  is  the  kernel  of  MAPLE®,  a 
system  developed  primarily  at  the  University  of  Waterloo,  Canada,  and,  more  recently,  at 
the  Eidgenossiche  Technische  Hochschule,  Zurich,  Switzerland.  MAPLE  is  marketed 
and  supported  by  Waterloo  Maple,  Inc. 

Maple  V  Release  4  was  incorporated  into  MATLAB®  5.  There  are  tow  toolboxes. 
The  basic  Symbolic  Math  Toolbox  is  a  collection  of  more  than  one  hundred  MATLAB® 
functions  that  provide  access  to  the  Maple  kernel  using  a  syntax  and  style  that  is  a  natural 
extension  of  the  MATLAB®  language.  The  basic  toolbox  also  allows  you  to  access 
functions  in  MAPLE®’ s  linear  algebra  package.  The  Extended  Symbolic  Math  Toolbox 
augments  this  functionality  to  include  access  to  all  nongraphics  MAPLE®  packages, 
MAPLE®  programming  features,  and  user  defined  procedures.  With  both  toolboxes,  you 
can  write  your  own  M-files  to  access  MAPLE®  functions  and  the  MAPLE®  workspace 
[Ref.  13]. 

In  order  to  determine  it  you  have  the  Extended  Symbolic  Math  Toolbox  a  simple 
check  on  the  MATLAB®  command  line  will  work.  The  syntax  is  as  follows: 

»maple(’with’,’numtheory’) 

This  will  either: 

1)  list  a  set  of  libraries  that  are  loaded  or 

2)  produce  the  following  error: 

???  ’with ’(...)’ requires  extended  symbolic  toolbox. 

Error  in  ==>  Paul  P.:Applications  :MATLAB  4.2a:Toolbox:Symbolic 
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lnstaller:symbolic:maplemex.mex 

Error  in  ==>  Paul  P.: Applications  :MATLAB  4.2a:Toolbox:Symbolic 
lnstaller:symbolic:maple.m 

On  line  84  ==>  [result, status]  =  maplemex(statement); 

More  information  about  obtaining  the  Extended  Symbolic  Math  Toolbox  may  be 
obtained  from  the  Math  works  website  f www.Mathworks.coml. 

This  thesis  is  written  assuming  the  reader  possesses  a  basic  understanding  of 
MATLAB®.  Additional  help  with  the  Symbolic  Math  Toolbox  may  be  obtained  in 
several  ways.  General  information  about  symbolic  functions  may  be  obtained  by  typing, 
»  help  sym /function. 

Here,  function,  is  the  name  of  an  “overloaded”  MATLAB®  numeric  function. 
This  method  provides  symbolic-specific  implementation  of  the  function,  using  the  same 
function  name.  This  provides  some  consistency  in  the  logically  programming  process 
and  the  functions  used  in  both  MATLAB®  and  MAPLE®.  For  example  diff  is  a  function 
that  is  used  numerically  as  well  as  symbolically  by  MATLAB®. 

»  help  diff 

will  produce  a  different  result  that 
»  help  sym/ diff 

Help  for  the  numeric  version  informs  you  that  the  function  name  is  “overloaded”. 
Overloaded  methods 

»  help  char/diff.m 
»  help  sym/diff.m 

Help  may  also  be  obtained  on  MAPLE®  commands  from  the  MATLAB®  work 

space. 
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»  mhelp  diff 

This  returns  the  help  page  for  the  MAPLE®  diff  function. 

B.  IDENTIFYING  SYMBOLIC  VARIABLES  IN  MATLAB® 

Using  the  Symbolic  Math  processing  capabilities  of  MATLAB®  is  different  than 
using  MAPLE®.  Of  particular  difference  is  the  requirement  that  a  variable  be  identified 
as  a  symbolic  variable  before  it  used  in  MATLAB®.  If  you  have  ever  received  the  error 
???Undefined  function  or  variable  Y 

you  have  fallen  victim  to  this  peculiarity.  The  exception  to  this  rule  is  if  a  variable  is  the 
result  of  an  equation  consisting  of  previously  identified  symbolic  variable.  Then  the 
variable  is  automatically  assumed  to  by  symbolic. 

A  symbolic  variable  may  be  identified  in  at  least  two  methods.  To  designate  one 
variable  as  a  symbolic  expression  use 
u  =  sym(‘u’) 
or 

b  =  sym(‘beta’). 

Here  ‘u’  is  established  as  a  symbolic  variable  u  and  b  as  beta.  These  variables  are 
defaulted  to  be  with  respect  to  x.  A  variable  my  also  be  assigned  to  an  expression, 
f  =  sym(‘a*xA2  +  b*x  +  c’) 

but  in  this  form  the  Symbolic  Math  Toolbox  does  not  create  variables  corresponding  to 
the  terms  of  the  expression,  a,  b,  c,  x.  To  perform  symbolic  math  operations  (e.g., 
integration,  differentiation,  substitution,  etc.)  on  f,  each  variable  must  be  created 


26 


explicitly.  This  may  be  accomplished  in  the  form  above  or  by  typing, 
sym  a  b  c  x. 

Again  these  variables  are  defaulted  to  be  with  respect  to  x.  The  independent  variables  is 
chosen  with  the  idea  that  they  are  typically  lower-case  letters  found  at  the  end  of  the 
Latin  alphabet  (e.g.,  x,  y,  or  z).  Therefore  the  closest  letter  to  ‘x’  alphabetically  is  chosen 
as  the  default  symbolic  variable.  If  there  are  two  equally  close,  the  letter  later  in  the 
alphabet  is  chosen.  Establishing  a  variable  with  respect  to  another  variable  may  be 
accomplished  in  the  following  format, 
zet  =  sym(‘zet(t)’). 

The  variable  zet  is  now  defined  with  respect  to  ‘t\  This  is  of  particular  importance  when 
performing  integration  and  differentiation  as  in  this  thesis.  Further  explanations  of 
defining  symbolic  variables  in  MATLAB®  may  be  found  in  reference  13  symbolic 
toolbox.  v 

C.  SUBSTITUTIONS 

When  processing  expressions  using  the  Symbolic  Math  Toolbox,  the  MAPLE® 

kernel,  accessed  by  MATLAB®,  does  not  “a  priori”,  treat  the  expression  as  a  number. 

MAPLE®  assumes  that  the  symbolic  variables  as  “a  priori”  indeterminate.  That  is,  they 

are  purely  formal  variables  with  no  mathematical  properties.  Consequently,  when 

calculating  an  expression,  the  variables  do  not  automatically  assume  a  numerical  value, 

as  assigned.  For  example, 

»  syms  x  y  z 
» 1  =  x  +  y  *  z 
»  x  =  1; 
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»  y  =  3; 
»  z  =  5; 
»1 


produces  the  answer, 

1  =  xA2+  y  *  z 
However, 

1  =  subs(l) 
produces 
16 

In  many  cases  a  specification  of  the  desired  form  of  output  is  required.  If  a  numerical 
output  is  desired  for  the  expression  and  it  requires  processing  beyond  basic  arithmetic, 
the  ‘subs’  function  may  be  incorporated  with  the  ‘double’  command  to  produce  a  double 
precision  array.  Double  precision  array  is  the  general  MATLAB®  workspace  parameter. 
Inquiry  into  the  class  of  each  variable  may  be  accomplished  by  typing, 

»  whos. 

A  list  of  the  Name,  Size,  Bytes,  and  Class  of  all  variable  accessible  to  the  workspace  are 
listed. 


D.  MAPLE®  FUNCTION 

The  ‘maple  function’  allows  the  user  to  access  MAPLE®  functions  directly  from 

(r) 

the  MATLAB  workspace.  This  function  takes  sym  objects,  strings,  and  doubles  as 
inputs  and  returns  a  symbolic  object,  character  string,  or  double  corresponding  to  the 
class  of  the  input. 
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It  is  important  to  know  the  syntax  of  the  calling  sequence  used  by  MAPLE® 
function  to  be  used.  This  information  may  be  obtained  by  using  ‘mhelp.’  The  general 
format  for  using  the  MAPLE®  function  is, 

A ( i ,  j  )  =  maple (' coeff EOM(i) , ddDOFq(j )) , 
where  ‘coeff  is  the  maple  function  to  be  used.  In  this  case,  ‘coeff  sets  A(i,j)  equal  to  the 
ddDOF(j)  coefficients  of  the  expression  EOM(i).  In  some  instances  the  format  of  the 
MAPLE®  function  is  not  consistent  with  form  used  in  MAPLE®.  For  example  the 
MAPLE®  function  ‘union’  requires  the  sets  to  be  evaluated  to  be  on  either  side  of  the 
function  instead  of  in  parentheses  after  it: 

setl  :=  setA  union  setB 
instead  of, 

setl  =  union(setA,setB). 

These  forms  are  not  conducive  to  using  the  maple  function  however  the  proper 
form  by  be  accomplished  by  using  strings  and  converting  the  syntax.  An  example  of  this 
is  provided  in  Appendix  I  with  the  MATLAB®  function  munion. 
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IY.  MAPLE®  TO  MATLAB®  CONVERSION 


A  brief  explanation  of  the  programming  tools  required  to  produce  the  results  of 
this  thesis  is  provided  above.  Unfortunately,  excessive  trial  and  error  went  into 
establishing  the  differences  in  the  programming  format  required  to  produce  the  desired 
output.  The  similarities  in  the  two  languages  proved  to  be  a  hindrance  rather  than 
beneficial. 

A.  DEFINING  THE  VARIABLES 

In  order  to  accomplish  the  derivation  of  the  equation  of  motion  using  Lagrange’s 
equation  in  MATLAB®  the  following  hurdles  need  to  be  crossed.  MAPLE®  allows 
variables  to  be  subindexed  i.e. 

.  yr  =  OLt  +  (f>k 

where  (j)k  is  subindexed  with  k,  as  is  commonly  used  in  mathematical  text  books.  In  this 

case  k  represent  an  individual  rotor  blade  and  for  the  three  bladed  model  k  =  three.  This 
mathematical  formulation  is  represented  as 
psi : =Omega*  t+Phi  [k] ; 

in  MAPLE®.  MATLAB®  uses  square  brackets,  [  ],  to  represent  an  array  or  matrix.  As  a 
result  there  is  no  way  to  represent  k  in  the  MATLAB®  workspace  as  it  is  used  in 
MAPLE  .  A  multidimensional  array  was  used  to  resolve  this  conflict.  A  third  index  is 
used  to  represent  the  values  of  each  blade.  In  this  form  MATLAB®  calculates  ‘psi’  three 
time  for  a  three  bladed  model.  The  syntax  in  this  form  is. 

Phi  =  [Phil  Phi2  Phi3]; 
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Psi  =  Omega*t  +  Phi; 

But  when  the  variable  is  called,  it  is  done  so  with  a  for  loop  from  1  to  k  (k  represents  the 
number  of  blades). 

for  i  =  l:k 

ml(:,:,i)  =  [cos(psi(i)  sin(psi(i))  0 
-sin(psi(i))  cos(psi(i))  0 
001]; 

end; 

ml  is  generated  3  times  with  a  different  psi  each  time.  This  results  is  psi  with  Phi,  Phi2, 
and  Phi3  in  each  successive  generation.  This  accomplishes  the  necessary  ‘place  holding’ 
required  for  each  blade’s  energy  terms  to  accomplish  the  derivation. 


B.  SET  MANIPULATION 

As  a  result  of  MAPLE® ’s  dependence  on  taking  derivatives  with  respect  to 
independent  variables,  the  time  function  notation,  which  represents  the  degrees  of 
freedom  and  the  time  rates  of  change  of  the  degrees  of  freedom,  must  be  converted  to 
independent  variable  notation  as  noted  in  reference  4.  To  accomplish  this  a  substitution 
of  sets  was  used  to  switch  the  necessary  terms  before  and  after  the  differentiation.  In 
MATLAB®  this  may  be  accomplished  in  two  ways.  The  first  is  similar  to  the  form  used 
in  the  original  derivation. 

DOFF  =  [ul ,  u2  ] ; 

DOFB  =  [zetl,  zet2,  zet3] ; 

DOF  =  [DOFF,  DOFB]; 

dDOF  =  dif f (DOF, t) ; 

ddDOF  =  dif f (dDOF, t) ; 

syms  ql  q2  q3  q4  q5  dql  dq2  dq3  dq4  dq5  ddql  ddq2  ddq3 
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ddq4  ddq5 


DOFq  =  [  ql  q2  q3  q4  q5  ]  ; 
dDOFq  =  [  dql  dq2  dq3  dq4  dq5] ; 
ddDOFq  =  [  ddql  ddq2  ddq3  ddq4  ddq5 ] ; 
setl  =  [  ]  ; 
set2  =  [  ]  ; 

for  i=l:5  %  size  of  DOFq  vector 

setl  =  [setl  maple ( ' ' , DOF(i) ,DOFq(i) )  ]; 
setl  =  [setl  maple  ('  '  =  w  , dDOF(i)  ,  dDOFq (i)  )  ]; 
setl  =  [setl  maple  ( ' '  =  w,  ddDOF  (i)  ,  ddDOFq(i)  )  ]; 
set2  =  [set2  maple( ' '= w ,DOFq(i) ,DOF(i) )  ];  . 
set2  =  [set2  maple (''= w , dDOFq (i) , dDOF (i) )  ]; 
set2  =  [set2  maple  ('  '  =  w  ,  ddDOFq  (i)  ,  ddDOF  (i)  )  ]; 

end 

setl  =  maple ( ' convert ' , setl, ' set ' ) ; 
set2  =  maple (' convert set2 set ')  ; 

This  method  defines  the  terms  required  as  symbolic  variables.  It  defines  the  sets 
and  then  uses  a  maple  function  to  set  the  terms  equal  to  each  other  for  later  substitution. 

It  then  uses  another  maple  function  to  convert  the  sets  into  a  form  that  is  recognizable  by 
MATLAB®. 

The  second  method  uses  the  ‘subs’  command  local  to  the  MATLAB®  workspace. 

By  establishing  the  necessary  variables  into  an  array  of  ‘old’  and  ‘new’  terms  they  can  be 

switched  term  for  term  in  a  given  variable  (array). 

DOFF  =  [ul,u2] ; 

DOFB  =  [zetl,  zet2,  zet3]; 

DOF  =  [DOFF,  DOFB]; 

' dDOF  =  diff (DOF, t) ; 
ddDOF  =  diff (dDOF, t) ; 
tempDOF  =  [DOF  dDOF  ddDOF] ; 

syms  ql  q2  q3  q4  q5  dql  dq2  dq3  dq4  dq5  ddql  ddq2  ddq3 
ddq4  ddq5 


DOFq  =  [  ql  q2  q3  q4  q5]  ; 


dDOFq  =  [  dql  dq2  dq3  dq4  dq5] ; 
ddDOFq  =  [  ddql  ddq2  ddq3  ddq4  ddq5 ] ; 
tempDOFq  =  [DOFq  dDOFq  ddDOFq] ; 

Temp  =  subs  (T,  dDOF,  dDOFq)  ; 

T,  ul,  u2,  zetl,  zet2,  zet3  are  predefined  arrays  and  tempDOFq  is  used  in  a  later  portion 
of  the  derivation.  The  speed  required  to  run  each  method  varies  depending  on  the 
operation  to  be  accomplished.  For  a  simple  substitution,  the  MATLAB®  ‘subs’  function 
proved  to  be  the  quicker  method  but  when  larger  symbolic  terms  were  involved  the  first 
method  proved  to  be  the  quicker  solution.  Therefore  an  optimization  was  done  to 
determine  the  optimum  combination  of  the  two  methods  in  order  to  minimize  the 
computational  time  required  to  complete  the  simulation. 

C.  STATE  DERIVATIVE  CALCULATION 

Using  the  Symbolic  Math  Toolbox  in  MATLAB  allows  calculation  of  the  state 
derivatives  without  first  creating  a  ‘B’  array,  as  is  necessary  with  the  MAPLE®  version. 
The  B  array  is  used  to  convert  the  output  to  C  or  Fortran  code,  before  it  is  placed  in  a 
MATLAB®  function.  By  avoiding  this  step,  both  the  MATLAB®  function  and  the 
SIMULINK®  function  can  call  the  state  derivative  directly  from  the  workspace.  This  is 
extremely  beneficial  in  that  it  eliminates  any  errors  associated  with  converting  generated 
code  to  a  format  that  usable  in  SIMULINK®..  It  eliminates  extra  steps  in  the  setup 
process,  helps  facilitate  simulations,  and  reduces  the  time  taken  to  convert  the  generated 
code. 
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D.  VARIABLE  ROTOR  BLADE  MODELING 


After  the  code  was  found  to  run  successfully  in  SIMULINK®,  it  was  further 
adjusted  to  allow  for  a  variation  in  the  number  of  rotor  blades.  In  the  MAPLE®  version 
this  was  self-contained  by  using  the  subindex  ‘k’.  In  the  MATLAB®  version,  without  the 
capability  of  subindexing,  each  array  was  built  up  as  a  temporary  array,  capable  of 
deriving  a  seven  bladed  model,  then  only  the  necessary  terms  were  used  to  complete  the 
derivation. 

E.  COMPLEX  MODEL 

The  complex  model  as  derived  by  Robinson,  [Ref.  4],  was  converted  into  a 
MATLAB®  function,  Appendix  E.  This  model  is  very  applicable  to  future  helicopter 
design  due  to  the  added  degrees  of  freedom.  It  is  the  complex  model  and  variations  of  it 
that  will  allow  the  designer  to  examine  the  stability  characteristics  of  rotor  systems 
without  expensive  and  time  consuming  tests  in  a  wind  tunnel.  The  flexibility  in  alteing 
the  degrees  of  freedom  and  accessing  the  equations  of  motion  with  each  time  step  are 
what  make  the  simulation  useful  beyond  the  calculations  performed  using  the  simple 
model.  Additionally  the  inclusion  of  nonlinear  terms  opens  a  whole  new  field  of 
dynamic  modeling  to  be  considered.  Some  of  the  possible  areas  of  exploration  are 
damper  springs,  duffing  springs,  and  the  stability  characteristics  of  varying  the  length 
rotor  blades. 

F.  MAPLE®  PROCEDURE 

Reference  13  describes  a  method  to  input  a  MAPLE®  symbolic  derivation  into 
the  MATLAB®  workspace.  Creating  a  source  file  using  a  ‘maple  procedure’  the 
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Extended  Symbolic  Math  Toolbox  is  capable  of  reading  the  source  file.  Unfortunately 
the  results  produced  little  to  no  gain  due  to  the  process  of  changing  the  source  code 
identifier.  This  process  was  further  complicated  due  to  windows  automatically 
associating  the  file  as  a  text  file  or  attaching  its  own  identifier  to  file.  Accessing  the 
sources  code  was  eventually  possible  by  changing  the  file  in  DOS  mode.  This  process 
could  possibly  be  simplified  using  a  UNIX  based  system  however  it  does  not  produce  the 
ultimately  desired  results  of  this  thesis. 
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V.  SIMULATION  MODEL 


A.  ORIGINAL  SIMULINK®  S-FUNCTION 

Introduction  of  the  rotor-fuselage  dynamics  into  the  Simulink®  simulation  is 

accomplished  through  the  use  of  a  Simulink®  S-function.  The  S-function  is  a  user 
defined  Simulink®  block  that  outputs  the  rotor  time  histories  as  functions  of  the  user 
defined  rotor  and  fuselage  parameters.  SIMULINK®  is  a  convenient  Graphical  User 
Interface  (GUI)  that  performs  simulation  time  integration,  allows  user  input  changes 
between  simulations,  and  provides  the  user  with  graphical  output. 

SIMULINK®  accesses  an  S-function  through  is  numerical  integration  routines. 
The  routines  make  calls  to  the  S-function  for  specific  information,  the  type  of  information 
returned  is  dependent  on  the  value  of  a  flag  variable  sent  by  the  integration  routine.  For 
example, 

flag  =  0  S-function  returns  sizes  of  parameters  and  initial  conditions, 
flag  =  1  S-function  returns  state  derivatives  dx/dt, 
flag  =  3  S-function  returns  outputs. 

The  section  of  the  S-function  which  computes  the  derivatives  at  each  time  step  is  the 
section  which  calls  the  equations  of  motion. 

The  Lagrangian  derivation  of  the  equations  of  motion  generated  the  equations  in 
the  form, 

F(x,x,x,t)  =  0 
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where  X  is  a  vector  of  displacement  degrees  of  freedom  of  the  system.  The  output  of  the 
derivation  was  further  manipulated  into  the  equivalent  form, 

[a(x,x,o]*  =  /(*,*,  0 

where  A  is  an  NxN  matrix  and  f  is  an  Nxl  vector,  with  N  =  number  of  degrees  of 
freedom  of  the  system.  This  is  possible  since  the  equations  are  quasi-linear  in  the  second 

derivative  terms,  i.e.,  no  terms  of  types  such  as  x ,  or  sin(  x ),  etc.  This  form  can  then  be 
transformed  from  N  second  order  equations  to  2N  first  order  equations  as  follows, 

X  =  vv 

*=M-7 

These  equations  are  evaluated  at  each  time  step  in  a  numerical  simulation  to  give 
the  state  derivatives.  This  study  uses  a  Runge-Kutta  algorithm  in  SEMULINK®  to  solve 
the  numerical  ordinary  differential  equations.  The  Runge-Kutta  algorithm  generally 
outperforms  other  schemes  for  systems  of  nonlinear  ordinary  differential  equations  which 
are  not  too  stiff  and  they  handle  discontinuities  well  [Ref.  14]. 

A  Coleman-like  3-bladed  hub-rotor  model  was  used  for  the  initial  simulation 
model.  It  provided  for  linear  stiffness  and  damping,  quadratic  damping,  and  cubic 
stiffness  in  the  blade  degrees  of  freedom.  Other,  simulations  for  damperless  rotor 
analysis  have  been  created  for  four-blade  and  five-blade  hub-rotor  simulations.  Table  2 
provides  a  break  down  of  the  variables  used  in  the  S-function  and  Figure  7  shows  the 
SIMULINK®  model  utilized  for  the  simple  model  simulation. 
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Rotor  Blade  Mass 

mb(1)  mb(2)  mb(3) 

mass 

Fuselage  effective  in  x  and  y  direction 

M(1)  M(2) 

mass 

Distance  from  hinge  to  center  of  mass  of  blade 

R 

length 

Rotor  Speed 

Omega 

rad/sec 

Hinge  Offset 

el 

length 

Angle  at  which  lead-lag  stops  engage 

z 

radians 

Azimuth  phase  angle  of  rotor  blade 

Phi(1)  Phi(2)  Phi (3) 

radians 

Lead-lag  linear  damping  coefficient 

Czeta(t)  Czeta(2)  Czeta(3)) 

moment/(rad/sec) 

Lead-lag  nonlinear  damping  coefficient 

Vzeta(l)  Vzeta(2)  Vzeta(3) 

Fuselage  linear  damping  coefficient  in  x  and  y  direction 

C(1)  c(2) 

force/(length/sec) 

Fuselage  nonlinear  damping  coefficient  in  x  and  y 
direction 

v(1)  v(2) 

Lead-lag  linear  spring  coefficient 

Ke(1)  Ke(2)  Ke(3) 

moment/rad 

Lead-lag  nonlinear  spring  coefficient 

Kd(1)  Kd(2)  Kd(3) 

Lead-lag  stop  spring  coefficient 

Ks(1)  Ks(2)  Ks(3) 

moment/rad 

Effective  fuselage  stiffness  in  the  x  and  y  directions 

K(1)  K(2) 

force/length 

Fuselage  states  initial  displacement  conditions 

xXi  xYi 

length 

Fuselage  states  initial  rate  conditions 

xrXi  xrYi 

length/sec 

Blade  states  initial  displacement  conditions 

xl  i  x2i  x3i 

rad 

Blade  states  initial  rate  conditions 

xrl  i  xr2i  xr3i  xr4i 

rad/sec 

Table  2.  Parameter  inputs  for  simple  model 


3- blade  rotor-pylon  (NL) 


Fuselage  Displacement 


Mux  Lead-lag  Displacement 


Figure  7.  SIMULINK®  diagram  of  NPS  simple  model 
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B.  MODIFICATIONS  TO  SIMULINK®  S-FUNCTION 

Of  significant  importance  in  this  and  future  versions  of  the  simulation  is 
MATLAB®’s  definition  of  variables.  From  the  original  MATLAB®  function  where  the  B 
matrix  was  input  as  optimized  C  or  Fortran  code,  the  variables  in  the  equations  of  motion 
were  in  array  format.  Using  the  Symbolic  Math  Toolbox  requires  the  variables  to  be 
defined  as  symbolic.  In  doing  so,  an  additional  substitution  step  is  required  to  solve  the 
state  derivatives  as  real  variables.  This  substitution  may  be  accomplished  in  either  the 
derivation  function  or  the  MATLAB®  S-function.  A  coordination  between  the  two 
functions  is  required  in  order  to  properly  derive  the  equations  of  motion  and  still  run  the 
simulation.  For  example  time  (t)  is  needed  to  derive  the  equations  of  motion,  however  it 
is  not  recognized  in  MATLAB®  unless  it  is  defined  as  a  symbolic  variable  prior  to  the 
derivation,  but  once  it  is  defined  as  a  symbolic  variable  an  actual  value  for  ‘t’  (each  time 
step)  may  not  be  substituted  in  to  the  simulation  model  until  it  is  redefined  as  a  ‘double 
array’  variable. 

Using  a  MATLAB®  function  to  derive  the  equations  of  motion  allows  the  model 
to  call  the  derivation  function  with  each  time  step.  Unfortunately  this  process  is  very 
time  consuming  and  severally  slows  the  simulation  process.  The  potential  exists  to  cut 
and  past  the  equations  of  motion  into  an  S-function,  in  order  to  reduce  the  simulation  run 
time.  This  method  was  not  used  in  an  attempt  to  maintain  the  integrity  of  the  function 
calling  process.  Future  thesis  may  have  the  opportunity  to  optimize  this  process.  It 
should  be  noted  that  the  cut  and  paste  method  would  still  significantly  reduce  the  time, 
labor,  and  possible  errors  associated  with  converting  the  MAPLE®  code. 
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As  an  alternative  method  to  deriving  the  equation  of  motion  with  each  time  step,  a 
data  file  was  created  with  the  state  derivatives  input  from  the  derivation  function.  The 
goal  of  this  setup  was  to  reduce  the  computational  time.  This  method  was  not  found  to  be 
as  effective  as  anticipated  because  the  process  still  required  the  derived  variables  to  be 
converted  to  ‘double  array’  variables  before  value  substitution  in  the  S-function 
workspace. 

C.  USING  THE  SIMULINK  MODEL 

Six  files  are  required  to  run  the  simulation.  For  the  three  bladed  simple  model  the 
files  are:  simple3.m,  helo3b5.m,  simple.m,  simple.mdl,  signum.m,  signuml.m,  abs.m. 
From  the  MATLAB®  command  window  type  in 
»  simple 

a  figure  will  appear  identical  to  Figure  7.  The  operator  may  begin  the  simulation  by 
“clicking”  on  the  arrow  (-►  )  on  the  menu  bar.  Two  figures  will  be  displayed.  Fuselage 
Displacement  and  Lead-lag  Displacement.  The  Fuselage  Displacement  figure  represents 
the  position  of  the  fuselage  in  the  x/y  plane  with  respect  time.  The  Lead-lag 
Displacement  figure  represents  each  blades  position  with  respect  to  time.  From  here  the 
operator  can  deduce  whether  the  model  is  stable  or  unstable. 

The  rotor  parameters  may  be  changed  by  “double-clicking”  on  the  box  titled  “3- 
blade  rotor-pylon  {NL}.  A  new  window  will  appear  containing  the  “Block  Parameters.” 
Each  of  these  terms  may  be  changed  to  effect  the  desired  response  from  the  model.  An 
explanation  of  these  terms  is  included  in  Table  2.  Further  explanation  of  using  the  NPS 
simulation  model  may  be  obtained  in  reference  17. 
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VI.  RESULTS 


A.  EQUATIONS  OF  MOTION 

The  equations  of  motion  calculated  for  the  simple  model  by  MATLAB®, 
Appendix  A,  are  equivalent  to  the  equations  calculated  by  Robinson  [Ref.  4],  Asa  result 
the  plots  produced  by  SIMULINK®  display  the  same  time  history  plots.  Figures  8,  9. 

The  derivation  of  the  equations  of  motion  using  MATLAB®  Symbolic  Math  Toolbox  are 
complete. 


Hub  Motion  in  Horizontal  Plane 


Hub  x-axis  displacement  (ft) 

Figure  8.  Fuselage  displacement  for  unstable  region 
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Blade  Lead-Lag  Angle, £  (deg) 


Lead-Lag  Displacement  Time  Histories 


0  0.1  0.2  0.3  0.4  0.5  0.6  0.7  0.8 


Time  (sec) 

Figure  9.  Lead-Lag  displacement  for  unstable  region 


VII.  CONCLUSIONS  AND  RECOMMENDATIONS 


Previous  work  by  Robinson  and  King  was  based  on  a  derivation  in  which  the 
equations  of  motion  were  derived  in  MAPLE®  then  converted  into  Fortran  or  C  where  the 
equations  of  motion  were  further  converted  in  a  MATLAB®  usable  format  (Method  A). 
The  focus  of  this  thesis  was  to  derive  the  equations  of  motion  using  the  Symbolic  Math 
Toolbox  in  MATLAB®(Method  B).  This  derivation  allows  direct  communication 
between  the  derivation  function  and  the  simulation  thus  eliminating  user  interface  in  the 
simulation  process. 

Based  on  the  simulation  process  utilized  in  this  thesis.  Method  B,  deriving  the 
equations  of  motion  using  MATLAB®  was  found  to  possess  both  advantages  and 
disadvantages. 

ADVANTAGES 

-  No  required  conversion  of  equations  of  motion  to  MATLAB®  readable  code 

-  Eliminates  “user”  interface 
Single  software  package 

DISADVANTAGES 

-  Time  consuming 

Symbolic  processing  format  is  not  as  usable  as  MAPLE® 

-  Does  not  provide  instant  feedback  to  equation  processing 

-  Requires  additional  substitution  step 
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TRADEOFF 


-  It  was  found,  using  Method  B,  that  eliminating  user  interface  from  the  simulation 
process  incurs  a  substantial  run  time  penalty 

-  With  minor  user  interface,  modification  to  Method  B  eliminates  the  time  penalty 
incurred 

-  The  minor  user  interface  still  reduces  the  potential  for  error  associated  with  code 
conversion 

-  Modification  to  Method  B  maintains  a  comparable  run  to  Method  A 

Expanding  on  the  ‘Tradeoff  Option”,  if  the  user  intends  to  utilize  the  NPS  modeler  in  the 
form  in  which  the  simulation  calls  the  derivation  function  internally.  Method  B,  then  a 
run  time  penalty  is  incurred.  However  if  modifications  to  this  process  are  acceptable,  the 
result  is  a  reduction  in  potential  errors  associated  with  code  conversion  and  a  savings  in 
the  time  required  to  perform  this  process. 

For  example,  potential  modification  to  the  simulation  process  used  in  Method  B  is 
the  option  of  “cutting”  the  resultant  matrices  from  the  MATLAB®  workspace  and 
“pasting”  them  into  the  simulation  workspace.  This  “cut  and  paste”  method  is  similar  to 
deriving  the  equations  of  motion  in  Method  A.  However,  it  eliminates  the  process  of 
converting  the  code  from  Fortran  or  C  to  MATLAB®  executable  code  thus  reducing 
“user”  interaction  and  eliminating  the  potential  for  user  contamination  of  data. 
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Future  emphasis  should  be  placed  on  optimizing  the  passing  of  variables  between 
the  MATLAB®  derivation  function  and  the  simulation  process.  Continued  analysis  of 
this  process  may  determine  a  more  effective  means  of  substituting  the  required  variables 
and  reduce  the  time  to  run  the  simulation.  With  further  optimization  of  the  NPS  modeler 
in  the  MATLAB®  format  and  the  continued  increase  in  computational  processing  speeds, 
eliminating  user  interaction  between  the  derivation  and  simulation  process  is  a  realistic 
goal. 
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APPENDIX  A.  MATLAB®  EQUATIONS  OF  MOTION  FOR  SIMPLE  INPLANE 
(COLEMAN)  MODEL  WITH  THREE  ROTOR  BLADES 


pretty(EOM) 

[  2 

[1/2  vl  dql  abs(l,  dql)  +  vl  dql  |  dql  |  +  Ml  ddql 

2 

-  2  mb2  R  %4  Omega  cos(q4)  dq4  -  mb2  R  %4  Omega  cos(q4) 

2  2 

-  mbl  R  %6  Omega  cos(q3)  -  mbl  R  %6  cos(q3)  dq3 

2 

-  mbl  R  %6  sin(q3)  ddq3  +  mb2  R  %3  Omega  sin(q4) 

2 

+  mb2  R  %3  sin(q4)  dq4  -  mb2  R  %3  cos(q4)  ddq4 
2  2 

-  mb3  R  %2  Omega  cos(q5)  -  mb3  R  %2  cos(q5)  dq5 

2 

-  mb3  R  %2  sin(q5)  ddq5  +  mbl  R  %5  Omega  sin(q3) 

2 

+  mbl  R  %5  sin(q3)  dq3  -  mb3  R  %l  cos(q5)  ddq5 

2 

-  mbl  R  %5  cos(q3)  ddq3  +  mb3  R  %1  Omega  sin(q5) 

2  2 

+  mb3  R  %1  sin(q5)  dq5  -  mb2  R  %4  cos(q4)  dq4 

-  mb2  R  %4  sin(q4)  ddq4  +  K1  ql  +  mbl  ddql  +  mb2  ddql  +  cl  dql 
+  mb3  ddql  +  2  mbl  R  %5  Omega  sin(q3)  dq3 

-  2  mbl  R  %6  Omega  cos(q3)  dq3  +  2  mb3  R  %1  Omega  sin(q5)  dq5 

-  2  mb3  R  %2  Omega  cos(q5)  dq5  +  2  mb2  R  %3  Omega  sin(q4)  dq4 
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2  2  2 

-  mb3  el  %2  Omega  -  mb2  el  %4  Omega  -  mbl  el  %6  Omega  , 

2  2 

1/2  v2  dq2  abs(l,  dq2)  +  v2  dq2  |  dq2  |  -  mb3  R  %l  Omega  cos(q5) 

2 

-  mbl  R  %6  sin(q3)  dq3  +  mbl  R  %6  cos(q3)  ddq3 

2 

-  mbl  R  %5  Omega  cos(q3)  -  mbl  R  %5  sin(q3)  ddq3 

2  2 

-  mb2  R  %3  Omega  cos(q4)  -  mbl  R  %5  cos(q3)  dq3 
+  mb2  R  %4  cos(q4)  ddq4  +  mb3  R  %2  cos(q5)  ddq5 

2  2 

-  mb3  R  %2  sin(q5)  dq5  -  mb3  R  %2  Omega  sin(q5) 

2 

-  mb3  R  %1  cos(q5)  dq5  -  mb2  R  %3  sin(q4)  ddq4 

2  2 

-  mbl  R  %6  Omega  sin(q3)  +  mbl  ddq2  -  mb2  R  %3  cos(q4)  dq4 

2 

-  mb3  R  %1  sin(q5)  ddq5  -  mb2  R  %4  sin(q4)  dq4 

2  2  2 

-  mb2  R  %4  Omega  sin(q4)  -  mbl  el  %5  Omega  -  mb3  el  %1  Omega 
+  mb3  ddq2  +  mb2  ddq2  -  2  mb3  R  %1  Omega  cos(q5)  dq5 

“  2  mbl  R  %5  Omega  cos(q3)  dq3  -  2  mbl  R  %6  Omega  sin(q3)  dq3 

-  2  mb3  R  %2  Omega  sin(q5)  dq5  -  2  mb2  R  %4  Omega  sin(q4)  dq4 

-  2  mb2  R  %3  Omega  cos(q4)  dq4  +  M2  ddq2  +  K2  q2  +  c2  dq2 

2  2 

-  mb2  el  %3  Omega  ,  mbl  el  Omega  R  sin(q3)  +  Czetal  dq3  +  Kel  q3 
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3  2  2 

+  Kdl  q3  +  Vzetal  dq3  abs(l,  dq3)  -  ul(t)  +  mbl  R  ddq3 

+  mbl  ddq2  R  %6  cos(q3)  -  mbl  ddq2  R  %5  sin(q3) 

-  mbl  ddql  R  %6  sin(q3)  -  mbl  ddql  R  %5  cos(q3) 

2 

+  2  Vzetal  dq3  |  dq3  | ,  Vzeta2  dq4  abs(l,  dq4) 

3 

+  2  Vzeta2  dq4  |  dq4  |  +  Kd2  q4  +  Ke2  q4  -  u2(t) 

2 

+  mb2  el  Omega  R  sin(q4)  -  mb2  ddql  R  %4  sin(q4) 

2 

-  mb2  ddql  R  %3  cos(q4)  +  mb2  R  ddq4  +  mb2  ddq2  R  %4  cos(q4) 

3 

-  mb2  ddq2  R  %3  sin(q4)  +  Czeta2  dq4  ,  Czeta3  dq5  +  Ke3  q5  +  Kd3  q5 

2 

-  u3  +  Vzeta3  dq5  abs(l,  dq5)  +  2  Vzeta3  dq5  |  dq5  | 

2 

-  mb3  ddql  R  %2  sin(q5)  -  mb3  ddql  R  %1  cos(q5)  +  mb3  R  ddq5 

+  mb3  ddq2  R  %2  cos(q5)  -  mb3  ddq2  R  %1  sin(q5) 

2  ] 

+  mb3  el  Omega  R  sin(q5)] 

%1  :=  sin(Omega  t  +  Phi3) 

%2  :=  cos(Omega  t  +  Phi3) 

%3  :=  sin(Omega  t  +  Phi2) 

%4  :=  cos(Omega  t  +  Phi2) 

%5  :=  sin  (Omega  t  +  Phil) 

%6  :=  cos(Omega  t  +  Phil) 
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APPENDIX  B.  MATLAB®  WORKSHEET  FOR  SIMPLE  INPLANE  (COLEMAN) 
MODEL  WITH  THREE  ROTOR  BLADES 


function  [Al,  fl]  =  simple3 ( ) 


Q.o.Q.9.9.Q.Q.o.e-9.@.9.o„9.9.9.9.e.9-e.ao.9.£.9„Q.a.9.o..g,9.9.9.9..g.9.9.£.9.g.9.9.£.9.9.s.9.9.e.g.9.9.e.9.9.g-9-9-9- 

o  o  "c  c  o  o  o  6  o  6  o  o  c  6  c  o  o  6  o  5  o  o  o  6  o  o  o  o  o  c  o  o  o  o  G  o  o  6  5  c  5  o  6  c  c  5  o  o  6  6  6  oo  ^  o  o  o  o  5 

%  EQUATIONS  OF  MOTION  FOR  A  HELICOPTER  IN  GROUND  RESONANCE 
%  CONSIDERING  ONLY  INPLANE  DEGREES  OF  FREEDOM 

aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaQ.aQ.Q.Q.aaQ.2,aaaQ.Q.Q.Q.aaQ.Q.aQ.aQ.aQ.Q. 

•b'6t'0'6'b-6'6'6t'6t'D'6t^tt'6'5'6'6-6t'6-6-6-6'6'5pO'S’6'6'6t-6’&'6t'6'8-6'6ttt^'6tt-6-6-6'5-6ttt 


clear 

clc 


%  Definition  of  variables  for  transformations 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms  Phil  Phi 2  Phi 3  Omega  t 


Phi  =  [Phil  Phi 2  Phi3]; 
psi  =  Omega*t  +  Phi; 


zetl  =  sym( ' zetl ( t) ' ) ; 
zet2  =  sym( ' zet2 ( t ) ' ) ; 
zet3  =  sym( ' zet3 (t) ' ) ; 
zet  =  [zetl  zet2  zet3]  ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  COORDINATE  TRANSFORMATIONS 

a  a  a  c>.  o„  a  a  a  a  a  a  a  a  9,  9. 9, 0.  o„  9, 9. 9.  9.  9.  9,  o.  o.  9  9  9.  9.  a  9.  9. 9,  9,  9. 9.  9  9.  9.  9.  9.  9  9  9,  999999 9. 9,  9, 9.  9. 9. 9, 
oooot'o'oott'ttt'ooottotttooooot'otoot'oottooott  o  o  o  oott)  totooottob 


for  i  =  1 :  3 

ml(:,:,i)  =  [cos(psi(i))  sin(psi(i))  0 
-sin(psi(i))  cos(psi(i))  0 
0  0  1]; 

m2(:,:,i)  =  [cos(zet(i))  sin(zet(i))  0 
-sin(zet(i))  cos(zet(i))  0 
0  0  1]; 

end 

ml  ( : )  =  ml ( : )  . ' ; 
m2  (  :  )  =  m2  (  : )  .  '  ; 

Q.Q.0.9,  ag.Q.0.0,  Q,Q.Q.g,Q.Q.g.Q.9.g.9.o.o.Q.g.9.Q.Q.Q.o.9.g.o.9.g.Q.9.Q.9.g.g.9.9.o.o,  9.Q.Q.ag.g.Q.9.9.9!,o,  g.9.9.9. 
-6-6-6'0-6-6-6-6'o-6-6-6-6'6-6-6-6-6-6'6-6-6-6'6-6-6-6-6-6-6-6-6-6-6-6-6-6-6-6'6-6-6-6-t>-6-6-6-6-6-6-6-6-6-6'o-6-6'o-6 

%  ENERGY  OF  ROTOR  BLADES 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Kinetic  energy  of  rotor  blade  (TBk) 
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syms  el  R  mbl  mb2  mb 3 


mb  =  [mbl  mb2  mb3 ] ; 
ul  =  sym( ' ul (t) ' ) ; 
u2  =  sym ( ' u2 ( t ) ' ) ; 

g  o,  g.  g  g  a  o,  o,  a  g,  o,  a  g  a  q,  g  g  g.  q,  q,  q,  a  a  o,  a  o,  q,  q,  a  q,  a  a  q,  c 
•g  o  'o  o  "o  o  ■o'o  o'o'o  o  o  o  "o  "o  *£>  g  g  g  g  g  o  o  g  g  g  g  g  g  g  g  g  g " 


9-  g  9,  a  g  o„  o,  a  o,  a  a  9.  a  a  o,  o.  g  a  a  o.  g  g  g  a. 
o  o  o  o  o  o  o  o  o  o  g  g  g  o  g  g  g  o  g  g  g  g  g  g 


rhoHI_I  =  [ul  u2  0]  ; 
rhoBuH  =  [el  0  0]  ; 
rhoPBd  =  [R  0  0]  ; 

for  i  =  1:3 

rhoBuH_I  =  rhoBuH*ml ( : , : , i) ; 

rhoPBd_I  =  rhoPBd*ml ( : , : , i ) *m2 ( : , : , i ) ; 

rho  =  simplify (rhoBuH_I  +  rhoPBd_I  +  rhoHI_I).'; 

V  (  : , : , i )  =  diff( rho , t ) . ' ; 

TBk ( : , : , i )  =  1/2 *mb ( i ) * (V( 1 , : , i ) ~2  +  V(2,:,i)A2); 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Potential  energy  or  rotor  blade  (UBk) 

g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  g  9,  g  g  g  o.  g  g  g  g  g  g  g  g  g  g  o,  g  g  g  g  g  o,  g  g  g  o,  9  a  g  g  a  a.  g  g  g  g  g.  a 
ggggggggg'6'6'6gggggggggg'og'6g'6g'6*6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6-6'6'5'6'6'g'g'g 


syms  Kel  Ke2  Ke3  Kdl  Kd2  Kd3  Ksl  Ks2  Ks3  z 
Ke  =  [Kel  Ke2  Ke3]; 

Kd  =  [Kdl  Kd2  Kd3 ] ; 


Ks  =  [Ksl  Ks2  Ks3 ] ; 

ggggggggggggggggggggggggggggg^gggggggggggggggg.ggggoi.ggQ,g.gggou  o, 

o  goo  go  g  o  6  6  6  6  g  g  o  6  o  0  5  o  6  5  o  6  go  'o  g  'o  6oq^>'o'o  t>  o'ooo  'o  '6  '6  "o  '6 15 15 1>  g  -6  g  g  '6  "6  g  g  '6 


for  i  =  1:3 

UBkl  =  1/2 *Ke (i) *zet (i) A2 ; 

UBk2  =  l/4*Kd(i) *zet (i) A4; 

UBk3  =  l/4*Ks  (i)  *signum(zet  (i) -z)  *  (zet  (i)  A2+zA2- 
2*zet  (i)  *z)  +l/4*Ks  (i)  *signum(zet  (i)+z)*(-zet(i)A2-zA2- 
2*zet (i) *z) +1/2 *Ks (i) *zet (i) A2  +  l/2*Ks (i ) *zA2 ; 

UBk(i)  =  UBkl  +  UBk2 ; 

end 


gggggggggggggggggggggggggggggggggggggggggggggggngggnnnngngg 

ggggggggggggggg'oggggggggggggggg'6'6gg'6'6'6'6'6'6'6'6'6'6-6'6'6'6’5'6-6'6-6g'6'6'6'6 

%  Dissapative  function  of  rotor  blades  (DBk) 

gggggggggggggggggggggggggggggggggggggggggggoggggggggggggggg 
'o'o'o'ogggggg'ooggog'gggggggg'o'o'6'6'6'6'6  6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6 

syms  DBl  DB2  DB3  Czetal  Czeta2  Czeta3  Vzetal  Vzeta2  Vzeta3 
Czeta  =  [Czetal  Czeta2  Czeta3]; 

Vzeta  =  [Vzetal  Vzeta2  Vzeta3]; 

ggggggggggggggggggggggggggggggggggggggggggggggggggggggggggg 

ggggggggggggggggggggggggggggggg'6g-^g'6'6'6'6'6-6'6'6'6'6'6'6'6'6'6-6gg^6'5'^'6g'5 


for  i  =  1 :  3 

DBk(i)  =  l/2*Czeta (i) * (diff (zet(i) ,t) )^2  + 
Vzeta  (i)  *  (diff  (zet  (i)  ,  t )  )  /v2*abs  (diff  (zet  (i)  ,  t)  )  ; 
end 


o,  9, 9^  o,  9. 9,  o,  9, 9, 9,  o,  o,  9. 9.  o,  o,  9^  c,  c,  q.  c,  o,  a  a  a  a  c,  a  o,  a  a  a  a  a  a  a  o,  a  a  a  o,  o,  g,  a  g.  o.  q.  o,  a  q.  a  o>  a  q,  a  a  o, 
o  o  o  Xs  o  0  o  o  o  o  o  o  0  o  t  t  o  ^  ^000  '0  0  '0000  ^  o  ^  *6  ”6  ^  ^  '6  'o  ^  ^  "6  ^  ^  ^  "6  *6  ^  'o  ti  "6  '6  'o  ^  'o  'o  "o  ^  *6 

%  ENERGY  OF  HUB 

9'^$'9'?-9'9'9'9'9'9-9-9-9'9'9'5'9'9'9'9-9-9-9'?^9'9'9'9'9-9'5>§-9-9-5.9-5.aQ.Q.aaaaaaao,o.ao,Q/o,ao,Q,ao, 
'o'o'6'b 'b'6'6'6  o  *6  t>  *o  *6  -6 "o  "6  -g  'o  '6 1>  'o  "6  “6  tS  "o  '6 -6 1>  'o  '6  'o  -*6 15  '6  "6  "5  '6  'o  15  'o  'o  'o  '6  "6  "5  "6  "6  "6  'o  'o  "5  '6  "6  'o  'o  'o  ^  'o 

%  Kinetic  energy  TH)  /  Potential  energy  (UH)  /  Dissapative 
function  (DH) 

^9-$'9'9'9'9'9'9-9'9'9'9-9'9'$'9'9'9-^9'9'9'9'9-9-Q'9'9-^9-9-Q-Q.9.a9-9-o,aaaaaaaQ.aaaao,ao.Q,o,Q.o.Q, 

t5^^l5l5l5^l5l5t5^^^t5t5l5t5l5l5l5l5l5^'o'6'6'6'6'6'6'6'£5'6'6'6''6'6'6'6'6'6'6'6'6'6l5'6'6't5'6'6#6'6'6'6'6,'6''6l5 

syms  K1  K2  cl  c2  vl  v2  Ml  M2 
c  =  [cl  c2 ]  ; 

v  =  [vl  v2 ]  ; 

M  =  [Ml  M2 ]  ; 

K  =  [Kl  K2]  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

TF  =  1/2*M(1) *dif f (ul , t ) A2  +  1/2*M  (2 ) *dif f (u2 , t) "2 ; 

UF  =  1/2  *K ( 1 ) *ul^2  +  1/2  *K (2 ) *u2^2 ; 

DF  =  l/2*c(l)*(diff (ul,t) )"2  + 

l/2*c  (2 )  *  (diff  (u2 ,  t )  )  /'2  +  l/2*v(l)  *  (diff (ul, t)  )  A2*abs  (diff  (ul 
,t) )+l/2*v(2) * (diff (u2,t) ) ~2*abs (diff (u2 , t) ) ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Generalized  forces  on  generalized  displacements 

9c  9r  9c  9c  9-  %  S'  9-  9-  9*  S'  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  S'  9-  9-  9-  9-  9-  9  9-  9-  9-  9-  9-  9-  9-  9-  9-  9~  9-  9-  9  9-  9-  ct  Q,  g^  g,  g.  g,  o.,  g,  g.  g,  g.  g.  g.  g,  g.  g, 

oouttyttt)atuototttottoto6uttu'ott'o'o'o  6  6  6  6'6'6'6'6'b'o'6'6'66t>'6'6t'6'6'6'6'6'6 

syms  FI  F2  F3  F4  F5  u3 
u  =  [ul  u2  u3 ] ; 

FI  =  0;  F2  =  0;  F3  =  ul;  F4  =  u2 ;  F5  =  u3 ; 

F  =  [FI  F2  F3  F4  F5]; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


DOFF  =  [ul , u2 ] ; 

DOFB  =  [zet (1) ,  zet (2 )  ,  zet(3)]; 
DOF  =  [DOFF,  DOFB] ; 
dDOF  =  diff (DOF, t); 
ddDOF  =  diff (dDOF, t) ; 


or  zet(l)  and  zet (2 


9c9c99999c999999c9999c9999999  9- 999999999  9- 999999  9- 9999999999aaaaaaa 
o  o  o  o  o  o  o  o  o  o  o  o'o'o  o  ot>t>t>t>t3t>t}'oi5t>t>t>i5t>^l5'g'6'o'6'6'6'6'6't>  6'6'6'6'6'6'6'6'6'6'6'’6'6'6'6'6'6'6 

%  DERIVATION  OF  EQUATIONS  OF  MOTION  USING  LAGRANGE'S 
EQUATION 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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syms  ql  q2  q3  q4  q5  dql  dq2  dq3  dq4  dq5  ddql  ddq2  ddq3  ddq4 
ddq5 

DOFq  =  [  ql  q2  q3  q4  q5  3 ; 
dDOFq  =  [  dql  dq2  dq3  dq4  dq5] ; 
ddDOFq  =  [  ddql  ddq2  ddq3  ddq4  ddq5] ; 


O'  <2*  2,  O'  O'  O'  O'  O'  O'  O,  O'  2,  ^  C,  0^0.  O,  O'  O,  O'  t 

'o'6'o'o'o'o'o'6'6yo'6'6'o'6'o'o'o'6'o'6'6'o'o'6i 


■9^^  9^  O' O'O'O'O'O'O.O'O'O'O'O'O' O' O'O'O'O'O'O'O'O'O'O'O'O.O' 

t'6'6'o;o^'6'6^'o'6'6'6'6'o'o:o'6'6'6'6'6'€'6'6;5'6'6''6'6^^i5;&'6 


9,9,^9,9,9,^9.9,9.9.o,9.9.9,9,o,o,9,9.9,9,9.9,9.9.aaaQ,o,o,Q.Q.aao,o,aaQ,Q,agQ.aaQ,Q,q,Q.aag,Q,Q,o.Q,a 
x>'d^^6^q^'oqoo^'6  6oo^^o^o^'o'o6  6'6  6'6'6'6i5'6  6'6'6'6'6'6'6'6'6'6''6'6'6;o'6'’6'6'6'6'6'6'6'6'6'6 

%  Substition  set  construction 


'o'o^'o'o'oo^'o^ooo'o'oto^o'o'o^bb'O'o'o'o'o'o^'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o^'o'o'o'o'o'o'o-'o'o'o 


setl  =  [ ]  ; 
set2  =  [  ]  ; 

for  i=l:5  %  size  of  DOFq  vector 

setl  =  [setl  maple ( ' ' ,DOF(i) ,DOFq(i) )  ]; 
setl  =  [setl  maple (' '=w ,dDOF(i) , dDOFq (i) )  ]; 
setl  =  [setl  maple  (''  =  w,  ddDOF  (i)  ,  ddDOFq  (i)  )  ]; 
set2  =  [set2  maple DOFq (i) , DOF (i ) )  ]; 
set2  =  [set2  maple  (  ,x  =  w,  dDOFq  (i)  ,  dDOF  (i)  )  ]; 
set2  =  [set2  maple  ddDOFq  (i)  ,  ddDOF  (i)  )  ]; 


end 

setl  =  maple (' convert setl set ') ; 
set2  =  maple (' convert ' , set2 set ') ; 

^tt^^^otootuoo'ootttt^ttt'oo'o'o'o'o'o'o'o'o'oo'o^'o'o'o'o'b'o'o^'o'o 


T  =  TF  ; 

U  =  UF; 
Dl  =  DF; 


[l,n] 
for  i 
T 
U 
Dl 

end 


=  size(DOFB); 

=  1:3 

T+TBk ( : ,  :  ,  i)  ; 
U+UBk ( i ) ; 

=  Dl+DBk ( i ) ; 


o,  O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  O' O' O' O' O' O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  0,  O'  O'  O'  9'  O'  C'  O'  O'  O'  O'  o,  o,  a  a  a  o,  q,  a  o,  o,  o,  a  q,  q.  o,  o,  q,  a  o,  q, 
'o  'o  'o  -o  'o  'o  *6  'o  'o  'o  *6  'o  *6  *6  'o  *o  *o  'o  'o  'o  "6  "o  'o  'o  x>  'o  'o  *o  *6  "o  'o  'o  'o  'o  'o  'o  'o  'o  'o  *6  'o  'o  *o  *6  -6  '6  *6  '6  *6  *6  *6  "6  "6  '6  *6  '6  *6  "6  '6 

%  Switch  of  time  dependent  terms 

9-  o,  o„  O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  O'  a  c.  a  p„  a  9-  9-  9,  o_  o„  o.  9,  a  a  q,  a  o,  a  a  q,  c.  a  o,  9,  q_  q,  a  q„  o.  c_  c„  oygag  q,  c„  o,o^o,g  q,  q,  a 
o  o  o  o  o  o  o  o  c  o  o  o  o  o  o  u  o  o  o  o  o  o  o  o  o  o  'o  *6  'b  'o  “t)  "o  "t)  "o  'o  “b  *6  'o  *6  'o  'o  ”6  'o  "o  "o  'o  "6  "o  'o  "o  ^  "o  "o  "o  "o  "6 

Temp  =  maple (' subs setl, T) ; 


ggggggggggggggggggggggggggggggggggggggggggggggggggggggggggg 

'6'6'6'6'6'6'6'o:5'6'6'6'6'6'6'6'6'6'6'6^'6'6'6'6'6'6'6'6'6'6'6'6'6'6-6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6^'6^'6'6'6^ 

%  Complete  EOM 

ggggggggggggggggggggggggggggggggggggggggggggggggggggggggggg 

'6'oo'o'o'6^'o'o'o^'o'6'6'o'o'6'o'o'o'o'6-o'o'6'o'o'o'6'6'o'6'6'c'6'6'6'6'6'6'6-'6'6'6'6'6'^'6'6'6'6'6''6'6'6'6'6"6'c 

[l,n]  =  size(DOF); 
for  i  =  1:5 
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tempi  =  diff (Temp, dDOFq(i) ) ; 
temp2  =  maple ( ' subs ' , set2 , tempi ) ; 
t emp3  =  di  f  f ( temp2 , t ) ; 

LI  =  maple (' subs ', setl, temp3 ) ; 

L2  =  diff (Temp, DOFq(i) ) ; 
tempL3  =  maple (' subs ', setl , U) ; 

L3  =  dif f (tempL3 , DOFq(i) ) ;  %  check  dimensions  of  DOFq 
tempL4  =  maple ( ' subs ' , setl , Dl ) ; 

L4  =  diff (tempL4, dDOFq(i) ) ;  %  check  dimension  of  DOFq 
EOM(i)  =  simplify (L1-L2+L3+L4-F (i) ) ;  %  check  dimension 

of  EOM  and  F 
end 


9.9,9,9.9.9,9.9,9,9.9,9,9.5,^9.9,^9,9,o,9,ao,ao,aaao,aaaaaaaaaaQ,Q,ao,aaaQ,aao.aaaaaaQ,a 

'6'd-6"6'o'6'6'o'6'6'6'6'6'6'D'6-^-6-6'o'6'6'6'6'6'6'6-6'6'6'6'6'6'6'6-6>6-6'6-6'6-6-6'6'6'6'6'6'6'6'6<6'6',6'6'6-6'6-6 


Elimination  of  excess  terms 


9.  9,  o.  o„  9, 9..  o„  9, 0.  9.  9, 9.  9. 9, 9. 9, 9, 9,  a  9.  9,  9,  9. 9. 9.  a  9.  a  9.  9.  9  a  9.  a  9. 9.  a  9-  a  a  a  a  a  a  a  a  a  q,  a  a  q,  a  o„  9  9.  a  a  a  a 
o  o  o  o  o  o  ^5  o  o  t>  o  o  "o  ^5  lo  o  15  6  "o  o  '0  6  '6  "6  -6 15  6  "6  o  "6  6  '6  6  'o  6  '6  6  'o  6  6  "6  '6  6  'o  o  '6  'o  'o  '6  i5  "0  ^  "6  'o  "6  'o  "6 


tempS  =  [  signuml (1, q3-z)  signuml ( 1 , q4-z )  signuml ( 1 , q5-z ) 

signuml (1, q3+z)  signuml ( 1 , q4+z )  signuml ( 1 , q5+z ) 

maple ( 'abs' , l,dq3)  maple (' abs ', 1, dq4)  maple ( 'abs' , l,dq5) 

maple (' abs ', 1, dql)  maple (' abs ', 1, dq2 )  ]; 

tempO  =  [00000000000]; 

sets  =  [  ]  ; 


for  i=l:ll 

sets  =  [sets  maple  (''  =  w,  tempS  (i)  ,  tempO  (i)  )  ]; 

end 

sets  =  maple (' convert ', sets ,' set ') ; 


9.9.9.9,99,9.999,9.9,9.9,9.9,o,9,9.o,9.ao.9aaaaao,o/9,o,o,o.o.Q,Q.9.aQ.aaaaao,o,D,g.g,ao,ao.aaaa 

x>^^^^^^t3'o'o'6'o'0'o^'o'o'o'o"o'o'6'o'o'o'o'6'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'o'oo'6'6'6'6'6'6'6'o'6'6'6'6'6 


GENERATE  A1  AND  fl  FROM  EOM 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


%  Determine  second  derivative  terms 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  i  =  l:n 
for  j  =  l:n 

A(i,j)  =  maple (' coeff ', EOM(i ), ddDOFq(j )) ; 

A(i,j)  =  maple (' subs ', sets, A (i , j )) ; 
end 

end 


Ax2dot  =  ddDOFq*A; 


9,9,9.9,5,9,9,9,9,^g,9.9,9,9,9,g,9,9,^5.a9,9,9,$,aao,$,9,ao,a9,aao,o,a6,aao,o,aaaaaao,o,ao,ao,ac, 

'o^'0'6o'6'o'6'o'5'o'&'6'6'6'o'6'6'o-t5'o'6'6'o'6-6'6'o'o'o'6'o'o'o-6'6'6'0'6'5'o'o'6'6'6'6'0'o'6'o'6'6'6'6'6'6'6'6'6 

%  Eliminate  second  derivative  terms  from  EOM 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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for  i  =  l:n 

f (i )  =  ( EOM ( i ) -Ax2dot ( i ) ) ; 
f  ( i )  =  maple ( ' subs ' , sets ,  f  ( i )  )  ; 
f  (i )  =  - (simplify (f (i) ))  ; 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%  Generate  state  vector 

o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  75  o  o  *6  *6  0  6  '6  'o  '6  6  6  6  '6  6  "6  '6  'o  "6  'o  '6  'o  'o  "o  'o  "o  'o  'o  'o  'o  '6'6  'o  'o  "o  '6 

syms  xl  x2  x3  x4  x5  x6  x7  x8  x9  xlO 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


Xldot  =  [xl  x2  x3  x4  x5]  ; 

XI temp  =  [x6  x7  x8  x9  xlO] ; 
Xl  =  [Xltemp  Xldot] ; 
DOFqtemp  =  [DOFq  dDOFq] ; 


'6'6'6'6'6"6'6'6'6'6'6'6 


9r9rSr9'9'9rgr9:9;9'9'9'9'9'9'9'9'9'9'9'9'9-9^9-9-9-9-Q.9.Q.Q.aaaaaaaaaaaaQ.Q,aQ. 
o  o  o  o  o  o ^  o  ox)  o  o^  o  o^5  ot5'c5'6  6'<5'6'6'6'6'6'6'6'6'6'6'6'6'6^'6'6'6''5'6'6'6'6'6'6'o 


%  Convert  terms  to  state  vector  form 


9:§:9'9'9'9'9'9-9'9'9'?'&'9'9-9'5-9-S'^9-9-$-&-9.aaQ.aaaaaaaaao. 
o  o  o'o  o'o  o  ^5  6  6  6  6^  x>66o  -6>  ^>6  '6  '66  'o'o  '6'6'o'6'6  'o  6S '6  yox>  '6  "o  'o 


setX  =  []  ; 
for  i=l:10 

setX  =  [setX  maple ( ' ' = ' ' , DOFqtemp (i) , Xl (i) )  ]; 

end 

setX  =  maple ( ' convert ' , setX, ' set ' ) ; 

A1  =  maple (' subs setX, A) ; 
fl  =  maple (' subs setX, f) ; 


9r0:§r§r9:gr9r9rg'P'P:9-9'9'9'9'9-9'9'9'9'9'9'9'9'9-9'9'9'9'9^9-Q-9-Q-9^o-Q-Q-Q-a(XQ.Q,Q,Q,Q,aQ.Q,Q„Q.Q,Q.Q,acLQ.a. 

%  Eliminate  variable  't'  from  input  'u' 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms  ul  u2  u3 
ut  =  [ul  u2  u3] ; 

Al  =  subs (Al,u,ut) ; 
fl  =  subs(fl,u,ut) ; 
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APPENDIX  C.  SIMULINK®  S-FUNCTION  FOR  SIMPLE  INPLANE 
(COLEMAN)  MODEL  WITH  THREE  ROTOR  BLADES 


function  [sys,  xO]  =  helo3b5 ( t , x, u, f lag, II , 12 , 13 , 14 , 15 , 16 ) 


%  function  [sys,  xO]  = 

helo3bA ( t , x,  u, flag, II, 12, 13,14,15,16! 


% 

s- 

S-function 

.  arguments : 

o 

% 

t 

= 

time 

% 

X 

= 

state  vector 

% 

u 

input  vector 

% 

flag 

= 

switch  used  by 

( simulation) 


function 

Sr 


routine  to  access  certain  parts  of  the  s- 
S-function  input  parameters: 


% 

% 


Q, 

Sr 


% 


11 

12 

13 

14 

15 


16 


[mb(l)  , mb ( 2  )  , mb  (3)  , M ( 1 )  , M ( 2 )  ] 

[R, Omega , el , z] 

[ Phi ( 1 ) , Phi ( 2 ) , Phi ( 3 ) ] 

[c  (1)  ,  c (2 ) , v(l) , v (2 ) , 

Czeta(l) ,Czeta(2) ,Czeta(3) , 
Vzeta ( 1) , Vzeta ( 2 ) , Vzeta ( 3 ) ] 

[  Ke  ( 1 )  ,  Ke  ( 2  )  ,  Ke  ( 3  )  , 

Kd ( 1 ) , Kd ( 2 ) , Kd ( 3 )  , 

Ks ( 1 ) , Ks ( 2 ) , Ks ( 3 ) , 

K(l) , K (2 ) ] 

[xrXi , xrYi , xrli , xr2i , xr3i , 
xXi , xYi , xli , x2i , x3i ] 


%  S-function  to  represent  dynamics  of  3  bladed  coupled 
rotor- 

%  fuselage  model  which  considers  only  inplane  degrees  of 
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%  freedom,  i.e.,  x  and  y  translational  fuselage  degrees  of 
freedom 

%  and  lead-lag  rotor  blade  degrees  of  freedom. 

% 

%  Explaination  of  variables: 

% - 


%  mb 

-> 

mass  of  blade 

%  M 

-> 

effective  mass  of  fuselage 

%  R 

of  mass 

-> 

distance  from  lead- lag  hinge  to  blade  center 

%  el 

-> 

blade  hinge  offset 

%  Omega 

-> 

rotor  speed 

%  z 

-> 

angle  at  which  blade  hits  stops 

%  Phi 

-> 

blade  phase  angle  w.r.t.  azimuth  postion 

%  c 

~> 

fuselage  linear  damping 

%  V 

-> 

fuselage  hydraulic  damping 

%  Czeta 

-> 

blade  linear  damping 

%  Vzeta 

-> 

blade  hydraulic  damping 

%  K 

-> 

effective  stiffness  of  fuselage  (landing 

gear  stiffness) 

%  Ke 

-> 

blade  elastic  spring  constant 

%  Kd 

-> 

blade  duffing  spring  constant 

%  Ks 

-> 

blade  stop  effective  spring  constant 

%  xr_  _i 

-> 

initial  rate 

%  x _ i 

a 

-> 

initial  displacement 

aaaQ.ao.aaa 
"o  -6  o  '6  "o  "6  'o  'o 

9'9'5'9-9'5'^5'9-5'9'9'?-9'9'9'9'5-9'9'9'9'9'c?'9'9'9-9'9'9-9-9-9-^Q'9-9-C,aaao,Q.o,aao,aao, 
o  o  o  o  o  ^5  x>  o  o  x>  75  to  o  o  75  15  t5  o  o  ^5  "b  o  o  'o  Taoooo  75  0  0  0  0  '0  To  75  To  75  "6  'b  75  '6  "0  'b  '6  '6  75  "6  "0 

%  Define 

input 

parameters 

$.0,0,0,  o.  q.  o,  a  a 
*o  *o  'o  "o  "o  'o  "o  "o  "o 

9. 9. 9, 9, 9, 9. 
•oooo  "O  "C 

.  9, 0. 9, 9. 9, 9, 9. 9, 9, 9, 9, 0. 9, 9, 9, 9, 9, 9, 0, 9, 0,  0, 0, 0,  a  a.  a  9, 9, 9, 9,  o„  q.  a  o„  9,  q,  9.  o,  a  o,  9,  cu  a 
i't>-6'b-9'&''o'6'6'6'6'6'6-6'6'6'6'b''o'6'6''o'6-6''o'6'o'6'6'6'6'6''o-'o'6''o^''6'6-'6'6"6'6''6'6 

%mbl=Il (1) ;mb2=Il (2) ;mb3=Il (3) ;M1=I1 (4) ;M2=I1 (5) ; 

%R=I2 (1) ;0mega=l2 (2) ;el=l2  (3) ;z=l2  (4)  ; 

%Phil=I3 (1) ; Phi2=l3 (2) ;Phi3=l3 (3) ; 

%cl=I4 (1) ;c2=I4 (2) ;vl=l4 (3) ;v2=l4 (4) ; 

%Czetal  =  I4 (5) ;Czeta2  =  I4 (6) ; Czeta3  =  l4 ( 7 )  ; 

%Vzetal=I4 (8) ;Vzeta2=l4 (9) ; Vzeta3=l4 ( 10 ) ; 

%Kel=I5 (1) ;Ke2=I5 (2) ;Ke3=I5 (3) ; 

%Kdl=I5 (4)  ; Kd2  =  15 (5) ;Kd3  =  I5 (6)  ; 

%Ksl=I5 (7) ;Ks2=I5 (8) ;Ks3=I5 (9) ; 

%K1=I5 (10) ;K2  =  I5 (11)  ; 

xrXi=l6 (1) ;xrYi=l6 (2) ;xrli=l6 (3) ;xr2i=l6 (4) ;xr3i=l6 (5) ; 
xXi=!6 (6) ;xYi=!6 (7) ;xli=l6 (8) ;x2i=!6 (9) ;x3i=!6 (10) ; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%. S-function  flag  conditionals 
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if  flag  ==  0 


sys= [10, 0, 10, 3 , 0, 0]  ; 

x0= [xrXi , xrYi ,  xrli , xr2i , xr3i , xXi , xYi , xli , x2i , x3i] ; 
elseif  abs(flag)  ==  1 
%  Calculate  derivatives 
[Al,  fl]  =  simple3; 

xl  =  x(l);  x2  =  x(2);  x3  =  x(3);  x4  =  x(4);  x5  =  x(5);  x6 

x  ( 6 )  ;  x7  =  x  ( 7 )  ;  x8  =  x  ( 8 )  ;  x9  =  x  ( 9 )  ;  xl  0  =  x  ( 10 )  ; 

ul  =  u(l)  ;  u2  =  u(2)  ;  u3  =  u(3)  ; 

x  =  subs (x) ; 
u  =  subs (u) ; 
fl  =  subs(fl); 

Al  =  subs(Al); 

sys  =  zeros (1, 2*5) ; 
sys (1:5)  =  Al\fl. ' ; 
sys (6:10)  =  x(l:5) ; 


%  Output  states 
elseif  abs(flag)  ==  3 


-0 

o. 


psi=Omega*t ; 
psil=psi+Phi (1) ; 
psi2=psi+Phi (2) ; 
psi3=psi+Phi ( 3 ) ; 

xcl=el*cos (psil) +R*cos (psil+x(8) ) ; 
xc2=el*cos (psi2 ) +R*cos (psi2+x ( 9 ) ) ; 
xc3=el*cos (psi3 ) +R*cos (psi3+x (10) ) ; 
ycl=el*sir± (psil)  +R*sin(psil+x(8)  )  ; 
yc2=el*sin (psi2 ) +R*sin (psi2+x ( 9 ) ) ; 
yc3=el*sin (psi3 ) +R*sin (psi3+x ( 10 ) ) ; 
xc=  (mb(l)  *xcl+mb(2)  *xc2+mb(3)  *xc3)  /  (sum(inb)  )  ; 
yc= (mb ( 1 ) *ycl+mb (2 ) *yc2+mb ( 3 ) *yc3 ) / ( sum (mb) ) ; 
Mag=sqrt (xcA2+ycA2 ) ; 
if  (xc<0&yc<0)  |  (yc  <  0) 

theta=pi+atan (yc/xc) ; 


o\°  o\°  o\° 


%  else 

%  theta=atan (yc/xc )  ; 

%  end 

%  gamma=theta+pi ; 

%  alphal= (psil/ (2*pi) -f loor (psil/ (2*pi ) ) ) *2*pi ; 

%  alpha2= (psi2/ (2*pi) -floor (psi2/ (2*pi) ) ) *2*pi ; 

%  alpha3= (psi3/ (2*pi) -floor (psi3/ (2*pi) ) ) *2*pi ; 

sys  (1 : 10)  =x; 

sys ( 11 ) =R*Mag*sin (gamma-alphal ) ; 
sys (12 ) =R*Mag*sin (gamma -alpha 2 ) ; 
sys (13 ) =R*Mag*sin (gamma -alpha 3 ) ; 
sys (14 ) = theta ; 

else 

sys  =  [ ] ; 

end 
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APPENDIX  D.  MATLAB®  WORKSHEET  FOR  SIMPLE  INPLANE  MOTION 
WITH  VARIABLE  NUMBERS  OF  BLADES 

function  [Al,  fl]  =  simple () 

99  99999-9  9  99'9999999999'9999999999999999Q,a.Q,Q.aaaaQ.Q,Q,Q,Q,Q,Q.Q, 

o  o  ot-b'o'o-o  <Dt'6'o'6'o-o'o'o^t'6'o'o^'6'6^^-6'o'6-6'o'6'6,c'6'6'-6'6'6'6'6^'6'6'6'6^'‘o'6'6'6l)'6'otio?1o 
%  EQUATIONS  OF  MOTION  FOR  A  HELICOPTER  IN  GROUND  RESONANCE 
%  CONSIDERING  ONLY  INPLANE  DEGREES  OF  FREEDOM 

9999999999999  9- 9999999999999999999999999  a  aaaaaaaaaaaaao,  cl  o.  a  a.  o. 

'6'o'6'o'o"o'o'6-6-6'o^>'o'6'6'6'6'6'6^'6'6'6'6'6^6'6'6'6'6-6-6-6-6'6-6'6-6'6'6'6%'6'6'6'6'6'6%%%%%%io%%%% 


%  mb 
%  M 
%  R 

of  mass 
%  el 
%  Omega 
%  z 
%  Phi 
%  c 

%  V 

%  Czeta 
%  Vzeta 
%  K 


-> 

-> 

-> 

-> 

-> 

-> 

-> 

-> 

-> 

-> 

-> 

-> 


mass  of  blade 

effective  mass  of  fuselage 

distance  from  lead- lag  hinge  to  blade  center 

blade  hinge  offset 
rotor  speed 

angle  at  which  blade  hits  stops 

blade  phase  angle  w.r.t.  azimuth  postion 

fuselage  linear  damping 

fuselage  hydraulic  damping 

blade  linear  damping 

blade  hydraulic  damping 

effective  stiffness  of  fuselage  (landing 


gear 

stiffness 

) 

%  Ke 

-> 

blade 

elastic 

spring 

constant 

%  Kd 

-> 

blade 

duffing 

spring 

constant 

%  Ks 
& 

-> 

blade 

stop  effective 

spring  constant 

"O 

a  a  a  a  a 
u  b  t  ho 

9  9-  9  9  9  9-  9-  9  9 
o  o  b  o  o  o  o  o  b 

9  9  9  9  9  9* 

O  O'O  'O  O  O' 

?„  o,  9, 9.  0.  0,  9.  Q.  a 
000  b>  b>  o  o  x>  "o 

a  a  a  a  a  a  a 
tt  o  o'o'o'o 

9'99'9999"999'9'99S'Q. 

b'obbbb  bbbbbbbbb 

clear 

clc 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Define  transformation  variables 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
k  =  3;  %  number  of  blades 


Phit  =  (Phil  Phi2  Phi3 
Phi  =  Phit(l:k); 
psi  =  Omega*t  +  Phi; 


Phi 4  Phi 5  Phi6  Phi.7]  ; 


%  blade  position  angle 


zetl  =  sym( '  zetl  ( t)  '  )  ;  zet2  =  syxn  ( '  zet2  ( t )  ' )  ;  zet3  = 
sym( ' zet3 ( t) ' ) ; 
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zet4  =  sym( ' zet4 ( t) ' ) ;  zet5  =  sym ( ' zet5 ( t ) ' ) ;  zet6  = 

sym( ' zet6 ( t ) ' ) ; 

zet7  =  sym{ ' zet7 ( t) ' ) ; 

zett  =  [zetl  zet2  zet3  zet4  zet5  zet6  zet7] ; 
zet  =  zett (1 :k) ; 


9-9'9'9'9-9'9'9'9'9'?'9'9'Q'9-Q'^^Q'5'Q'9'C,Q,Q,ao,Q.o,o/o,Q,o<o,Q.Q,o,aQ,aaaQ, 
o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  'o  'o  'o  “6  "o  'o  "o  "o  "o  o  'o  'o  'o  "o  'o  'o  'b  "o  "o  "o  *6  "6  *6  *6  *6  % 


%  COORDINATE  TRANSFORMATIONS 


%  %  %  Sc  %  Sc  9*  9-  9-  9-  9-  3-  9-  ^  9-  9-  Si-  9-  9-  9-  9-  9-  2-  9-  9-  9-  9*  9-  9-  a  9-  a  9-  a  a  a  a  cu  a  a  o,  q.  q.  q.  a  a  a  a  a  a  a  a  a  o,  a  q.  o„  q. 
o  o  o  o  o  o  o  o  o  o  o  o^>  o  o  o  o  o  o  o  o  o  too  o  o  t  o  o  6  6  75  o  o  o  "o  o  o  "6  ~o  "o  t>  "o  'o  'o'o'o'o  "6  '6  '6  75  'o  "6  '6  "6  "5  '6 


for  i  =  l:k 

ml(:,:,i)  =  [cos(psi(i))  sin(psi(i))  0 
-sin(psi(i))  cos(psi(i))  0 

001];  %  first  transformation 

m2(:,:,i)  =  [cos(zet(i))  sin(zet(i))  0 
-sin(zet(i))  cos(zet(i))  0 

001];  %  second  transformation 

end 

ml ( : )  =  ml 

m2  (  :  )  =  m2  ( :  )  .  '  ; 


9,9,^9,9,9,9,9,9,9,ao,aaaaQ,Q,aaao,aaQ.o,aQ,aQ.aQ,Q,Q,Q,o,Q.aao,Q,Q,Q,Q,oJ.Q,Q,Q,o.aQ,aoJ<o,o  o  0,0  o 

Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  ""Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  Q  "Q  Q  'Q  ^  "q  'q  'q  ^  ^  ^  'q  ^  ^  ^  ^  ^  ^  ^ 

%  ENERGY  OF  ROTOR  BLADES 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Kinetic  energy  of  rotor  blade  (TBk) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms  el  R  mbl  mb2  mb 3  mb4  mb5  mb6  mb7 
mbt  =  [mbl  mb2  mb3  mb4.  mb5  mb6  mb7]  ; 
mb  =  mbt (1 : k) ; 


ul  =  sym('uKt)');  u2  =  sym(  'u2  ( t)  '  )  ;  u3  =  sym( '  u3  ( t )  '  )  ;  u4 
=  sym ( ' u4 ( t ) ' ) ; 


u5  =  sym('u5(t) 


u6  =  sym( 'u6 (t) ' ) ;  u7  = 


^,9,9.9.^5.^a9,aaaaaaaaaaaaQ-aaaaaaaaag.naac 

'6'6'6-6'6'6-6^'6'6'6-6'6'6'6'6-6-6'6^6'6'6-6'6'6-6'6^6'6'6-6'6'6'6'6'i 


sym  (  '  u7  ( t )  '  )  ; 

%%%%%%%%%%%%%%%% 


rhoHI_I  =  [ul  u2  0]  ; 
rhoBuH  =  [el  0  0]  ; 
rhoPBd  =  [R  0  0] ; 


for  i  =  l:k 

rhoBuH_I  =  rhoBuH*ml ( : , : , i) ; 

rhoPBd_I  =  rhoPBd*ml ( : , :  ,  i ) *m2 ( : , : , i ) ; 

rho  =  simplify (rhoBuH_I  +  rhoPBd_I  +  rhoHI_I) . ' ; 

V ( : , : , i )  =  diff( rho , t ) . ' ; 

TBk ( : , : , i )  =  l/2*mb(i) * (V(l, : , i) A2  +  V(2,:,i)A2); 

end 
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o.  o..  9.  ct  q.  o-  9'  9-  9-  9-  9-  9-  9-  9-  o.  9. 9.  9  cu  9,  a  a  a  a  q,  9-  9-  a  Q,  a  o,  q,  o.  o,  o,  o.  o.  a  a  o,  g,  a  o,  n  a  o,  o.  a  o.  a  q,  o.  o.  o.  a  ou 

'6'6'6i'6'o'6'6^'6'6'6'6w6'6'6'6'6'6'6'6'6'6'6'6-6'6'6'6'6"o'6-6'6'b'6'6'6'5'6'6',6'b'6'6'6'6'6'fc>^6'6'6'6'b'6'6'6<6'6'b 

%  Potential  energy  or  rotor  blade  (UBk) 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms  Kel  Ke2  Ke3  Ke4  Ke5  Ke6  Ke7 

syms  Kdl  Kd2  Kd3  Kd4  Kd5  Kd6  Kd7 

syms  Ksl  Ks2  Ks3  Ks4  Ks5  Ks6  Ks7  z 

Ket  =  [Kel  Ke2  Ke3  Ke4  Ke5  Ke6  Ke7]; 


Ke  =  Ket (l:k) ; 

Kdt  =  [Kdl  Kd2  Kd3  Kd4  Kd5  Kd6  Kd7] ; 

Kd  =  Kdt (l:k) ; 

Kst  =  [Ksl  Ks2  Ks3  Ks4  Ks5  Ks6  Ks7]; 

Ks  =  Kst (l:k) ; 

aaaQ.aaaa9,aaaa9,9.aaaa9,aQ.aao.aaaQ,aaaQ.aaaaQ,Q.Q.aaaQ.Q,Q.aaag.aaaQ,aaaaQ, 

'o<of'ooboo'obbtt'o'6'6'D,6'obotbb’6'6'o'o't>'6'6'oo'o'o'6'6'o'bt)'o'6'o'ot)'o'b'o'6'6‘o,6'6''o'6'5'6'6'6 


for  i  =  l:k 

UBkl  =  l/2*Ke(i) *zet(i)^2; 

UBk2  =  l/4*Kd(i) *zet (i) ^4; 

%UBk3  =  l/4*Ks  (  i )  *signum(  zet  { i ) -z)  *  (  zet  ( i )  ~2  +  z/N2- 
2  *zet ( i ) *z)+l/4*Ks(i) *signum ( zet ( i ) +z ) * ( -zet (i) /'2-z~2- 
2*zet(i)*z)+l/2*Ks(i)*zet(i)~2  +  l/2*Ks (i) *z"2 ; 

UBk(i)  =  UBkl  +  UBk2 ; 

end 


Q,aaaao,Q,Q,Q,aaaaaQ,aap/aQ,aQ,Q.aaaaaQ,aaQ,Q.aQ,aQ,aQ,g.o,aQ,Q,aQ.Q,aQ,oyQ,Q,o,aQ,o,ao,a 

■o'oxto'ooootot'o'oo’ooo'ooooot-o'oo'o'ct'o'o'o'o'o't/o'o,t'6'o'o'o'o'6'o''o'ott'6'ot)'6'6''ot''o'6 


%  Dissapative  function  of  rotor  blades  (DBk) 


syms  DBl  DB2  DB3  DB4  DB5  DB6  DB7 

syms  Czetal  Czeta2  Czeta3  Czeta4  Czeta5  Czeta6  Czeta7 
syms  Vzetal  Vzeta2  Vzeta3  Vzeta4  Vzeta5  Vzeta6  Vzeta7 
Czetat  =  [Czetal  Czeta2  Czeta3  Czeta4  Czeta5  Czeta6 
Czeta7] ; 

Czeta  =  Czetat(l:k); 

Vzetat  =  [Vzetal  Vzeta2  Vzeta3  Vzeta4  Vzeta5  Vzeta6 
Vzeta7 ] ; 

Vzeta  =  Vzetat(l:k); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


for  i  =  1:3 

DBk(i)  =  l/2*Czeta (i) * (dif f (zet (i) , t) ) A2  + 
Vzeta  (i)  *  (dif  f  (zet  (i)  ,  t)  )  ~2*abs  (dif  f  (zet  (i)  ,  t)  )  ; 
end 
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%  ENERGY  OF  HUB 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Kinetic  energy  TH)  /  Potential  energy  (UH)  /  Dissapative 
function  (DH) 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms 

K1 

K2  cl 

c  = 

[Cl 

c2]  ; 

v  = 

[Vl 

v2  ]  ; 

M  = 

[Ml 

M2]  ; 

K  = 

[K1 

K2]  ; 

%%%% 

%%%% 

a  Q.  <9  Ou  Q,  Q, 
-6  "o  '6  'o  "o 

TF  = 

1/2*M ( 1) * 

UF  = 

1/2 

*K(1) * 

DF  = 

1/2 

*c (1) * 

1/2* 

c  (2 ) 

* (diff 

,  t)  ) 

+  1/2 

*v(2)  * 

9-  9"  9"  9-  9-  9"  9"  9'  9-  9--  9-  °-  Q*  9-  9-  9-  9-  9*  Q,  5-  9  9  Q.  9  9  9  o„  o„  o,  o_  o,  cx  q.  o„  o,  o,  o>  o_  o  o,  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  n  n 
o  o  o  o  o  o  o  o  o  o  o'o'o  o  o  o  o  o  "o  o  o  'o 'o  'o  'o  'o  'o  'o 'o  'o  ~b  'o  "o  'o  'o  'o  'o  'o  '6  ~o  'o  'o  'o  ~o 'o  '6  'o^  "6 'o  '6  'o  'o%  'S  ~6  'o  'o 


%%%%%<k§c%9;2'2'9'9''9'99999  9  9  9-'99  9  9  99999  9  9  99  9  9999  9  999<9  99Q,Q.o,.o  o  i 
o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o'b'o'o'o'0'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6 

%  Generalized  forces  on  generalized  displacements 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

syms  FI  F2  F3  F4  F5  u3 

ut  =  [ul  u2  u3  u4  u5  u6  u7] ; 

u  =  ut (1 :k) ; 

FI  =  0;  F2  =  0;  F3  =  ul;  F4  =  u2 ;  F5  =  u3 ;  F6  =  u4;  F7  = 
u5 ;  F8  =  u6 ;  F9  =  u7 ; 

F  =  [FI  F2  u]  ; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

DOFF  =  [ul , u2 ] ; 

DOFB  =  zet;  %  or  zet(l)  and  zet(2) 

DOF  =  [DOFF,  DOFB]; 
dDOF  =  di f  f ( DOF ,  t )  ; 
ddDOF  =  diff (dDOF, t) ; 

[l,n]  =  size(DOF); 


o  o  o  o  o  o  o  o  o  o  o  o  o'o'o'o'o'o'o'o'o'o'6'6'5'6'6'6'6'6'6'6'6-6'6'6'6-6'6'6'6'6'6'6'6'6'6%'6'6'6'6"6%%%%%% 

%  DERIVATION  OF  EQUATIONS  OF  MOTION  USING  LAGRANGE'S 
EQUATION 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

syms  ql  q2  q3  q4  q5  q6  q7  q8  q9 


syms  dql  dq2  dq3  dq4  dq5  dq6  dq7  dq8  dq9 

syms  ddql  ddq2  ddq3  ddq4  ddq5  ddq6  ddq7  ddq8  ddq9 

DOFqt  =  [  ql  q2  q3  q4  q5  q6  q7  q8  q9]  ; 

DOFq  =  DOFqt (l:n); 

dDOFqt  =  [  dql  dq2  dq3  dq4  dq5  dq6  dq7  dq8  dq9] ; 
dDOFq  =  dDOFqt ( l:n); 

ddDOFqt  =  [  ddql  ddq2  ddq3  ddq4  ddq5  ddq6  ddq7  ddq8  ddq9 ] ; 
ddDOFq  =  ddDOFqt (l:n); 


%  Substution  set  construction 


setl  =  [ ] ; 
set2  =  [ ] ; 

for  i=l:n  %  size  of  DOFq  vector 

setl  =  [setl  maple ( ' '=' ' ,DOF(i) ,DOFq(i) )  ]; 
setl  =  [setl  maple  (''  =  w,  dDOF  (i)  ,  dDOFq(i)  )  ]; 
setl  =  [setl  maple  (''  =  w,  ddDOF  (i)  ,  ddDOFq  (i)  )  ]; 
set2  =  [set2  maple ( ' '=' ' ,DOFq(i) ,DOF(i) )  ]; 
set2  =  [set2  maple( ' '=w ,dDOFq(i) ,dDOF(i) )  ]; 
set2  =  [set2  maple (''= w , ddDOFq (i) , ddDOF ( i ) )  ]; 

end 

setl  =  maple ( ' convert ' , setl , ' set ' ) ; 
set2  =  maple ( ' convert ' , set2 , ' set ' ) ; 


'o''o''o'6'6'b'6'o'6''o''o'6''o'6'6'^'6'b'6'o'6'6'6't)'b'6 


T  =  TF; 
U  =  UF; 
D1  =  DF; 


for  i  =  l:k 

T  =  T+TBk ( : ,  : , i )  ; 
U  =  U+UBk ( i ) ; 

Dl  =  Dl+DBk ( i )  ; 

end 


o.aaaaaaaaaaaaa9,aaQ,aact9,aaaaaQ,aaaaQ.aaaaaao.ao,aQ,o,o,aaaaaao.aaaaQ, 

'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O^'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'OO'b'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 

%  Switch  of  time  dependent  terms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

Temp  =  maple (' subs setl , T) ; 


9'9'9'9'9'9'5'9'9'^'9'P''9-9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9'9-9'9'9'9'9'-9'9'9'9'9'0-'$'9'c>''c^9'0-'9'Q'9'9'9'9'9'9' 
b  ^ooo  b  b  b  b  b  h  'o  b  b  "o  "o  b  b  b  "o  b  b  'o  h  o  b  "6  6  •o  b  o  "6  -6  b  b  b  o  o  "6  -66b6  'o  “6  "o  b  o  "6  'o  '6  0  ^  "6  'o  -6  "6  -6 

%  Complete  EOM 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

[l,n]  =  size(DOF); 
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for  i  =  l:n 

tempi  =  dif f ( Temp , dDOFq ( i ) ) ; 
temp 2  =  maple ( ' subs ' , set2 , tempi ) ; 
temp3  =  dif f ( temp2 , t ) ; 

LI  =  maple ( ' subs ' , setl , temp3 ) ; 

L2  =  di f  f ( Temp , DOFq ( i ) ) ; 
tempL3  =  maple (' subs setl , U) ; 

L3  =  diff (tempL3 ,DOFq(i) ) ;  %  check  dimensions  of  DOFq 
tempL4  =  maple (' subs setl , Dl ) ; 

L4  =  dif f (tempL4, dDOFq(i) ) ;  %  check  dimension  of  DOFq 
EOM(i)  =  simplify (L1-L2+L3+L4-F ( i )) ;  %  check  dimension 

of  EOM  and  F 
end 


^  S:  S;  S;  9-  9-  9'  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  ^  9-  9-  9-  9-  9- a  °-  9-  9-  c,  a  a  o,  a  a  a  a  Q,  q,  o,  o,  q,  o,  a  o  cl  a, 

o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o'o'o'o'o'o'o'o  o  o  "6  ^  "6  t>  "o  "o  -6  "o  "o  "o  "o  'o  -6  "o  "6  'S  "6  "6  ^6  -6  "6  "6  "o  % 


%  Elimination  of  excess  baggage 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


tempSmin  =  [  signuml  ( 1 ,  q3-z )  signuml ( 1 , q4-z )  signuml ( 1 , q5- 
z)  signuml (l,q6-z)  signuml (l,q7-z)  signuml ( 1, q8-z) 
signuml (1, q9-z)  ]; 

tempSpls  =  [  signuml ( 1 , q3+z )  signuml ( 1 , q4+z ) 

signuml (1, q5+z)  signuml ( 1 , q6+z )  signuml ( 1 , q7+z ) 

signuml (1, q8+z)  signuml (1 , q9+z )  ]; 

tempSabs  =  [  maple (' abs ', 1 , dql )  maple (' abs 1 , dq2 ) 

maple ( 'abs' , l,dq3)  maple ( 'abs' , 1, dq4)  maple (' abs ', 1 , dq5 ) 

maple ( 'abs' , l,dq6)  maple (' abs ', 1 , dq7 )  maple (' abs ', 1 , dq8 ) 

maple ( ' abs ' , 1 , dq9 ) ] ; 

tempS  =  [tempSmin (1 :n)  tempSpls ( 1 :n)  tempSabs ( 1 : n)  ]; 

[g,p]  =  size(tempS); 
tempO  =  zeros(l,p); 
sets  =  [ ] ; 


for  i=l:p 

sets  =  [sets  maple (''= w, tempS ( i ), tempO ( i ) )  ]; 

end 

sets  =  maple (' convert ', sets, 'set'); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


GENERATE  Al  AND  fl  FROM  EOM 


9-9-9-9-9-9-9-9-9-9-9-9-9-9-9.9-9- 
o  o  o  o  o  o  o  o  o  o  o  6  o  o  6  6 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Determine  second  derivative  terms 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

for  i  =  l:n 


for  j  =  l:n 

A ( i , j )  =  maple ( 'coeff ' ,E0M(i) ,ddD0Fq(j) ) ; 
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A(i,j)  =  maple (' subs ', sets, A( i , j )) ; 
end 

end 


Ax2dot  =  ddDOFq*A; 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Eliminate  second  derivative  terms  from  EOM 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


for  i  =  l:n 

f(i)  =  (EOM(i) -Ax2dot (i) ) ; 
f(i)  =  maple (' subs ', sets, f (i) ) ; 
f(i)  =  - (simplify (f (i) )) ; 

end 


-6 -6 15 -6  'o  t>  o  o  "6  "6  '5  -6  'o  'o  "6  -6  "6  -6  15  'o  15  -"o  "6  '6  “6  ^6  'o  'o  'o  ^6  "o  'o  "6  -*6  -6  "o  'o  "6  ^  "6  'o  "6  "6  '6  "6  "6  '6  "6 15  %  %  t>  %  %  %  %  % 


Generate  state  vector 


o,9,9,9,o,ci.a^c,o,Q,ao,o.Q,Q,aanaao.ao.ao,a 

'6'6'6'o'6'6’€'6'o'o'6'6'6'6'6'o'6'6'6'6'6'6'6'6-6-6'6 


'6'6'6'6'6'6'6'6'6'6'6,b'6'6'6'6'6'6'6'6'6'6'6'6'6'6 


syms  xl  x2  x3  x4  x5  x6  x7  x8  x9  xlO  xll  xl2  xl3  xl4 

§:§;§:S:§cS:S:5'^9'9'9'9'9-9'9'9'9-3'9'9'9-9'9'9-9‘9'9'g'5'5-5-5'9'9'9'9<9<Q'3-9'9.5-9-9-aa9.9.ag.aaaagaaa 

'b'o'o'o'o-t>x>l3x>'b'b'o'o'o'o^'o'o'o'o'o^o'o-o'6'o'6'6'6'6'6"o'6'6''6'^'6"6'6'6'6'6'6'6'6'6'6^'6'6^-6'6'6'6'6'6'6^ 


Xldot  =  [xl  x2  x3  x4  x5] ; 
Xltemp  =  [x6  x7  x8  x9  xlO] ; 
Xl  =  [Xltemp  Xldot] ; 
DOFqtemp  =  [DOFq  dDOFq] ; 


9-  9-  ^  9-  9'  9-  9-  9-  9*  ^  9-  9-  9-  9^  Q-  5^  5-  5-  ^  ^  9-  9-  0-  o,  a  a  q,  a  a  o,  a  a  q.  o,  q,  o,  a  o.  a,  a  q.  o,  q,  o,  a,  a, 

o  -6  -6  -6  "o  'o  *6  "6  "o  'o  -6  ■$  S  'o  "o  "o  'o  "o  *6  "o  *b  t>  "6  "o  *6  *6  -6  -6  "6  "6  -6  -6  •'6  t>  "6  "6  *6  *6  '6  *6  -6  "o  -6  -6  "6  '6  '6  "6  *6  '6  '6  *6  '6  "6  %  %  %  %  % 

%  Convert  terms  to  state  vector  form 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

setX  =  [  ]  ; 
for  i=l:2*n 


setX  =  [setX  maple DOFqtemp (i ) , Xl (i ) )  ]; 

end 

setX  =  maple (' convert ', setX, 'set'); 


A1  =  maple (' subs ', setX, A) ; 
fl  =  maple (' subs ', setX, f); 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  Eliminate  variable  't'  from  input  'u' 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


syms  ul  u2  u3 
ut  =  [ul  u2  u3]; 
A1  =  subs(Al,u,ut); 
f  1  =  subs(f  1  ,u,ut); 
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APPENDIX  E.  MATLAB®  WORKSHEET  FOR  COMPLEX  (STRAUB)  MODEL 
WITH  VARIABLE  NUMBER  OF  ROTOR  BLADES 

function  [Al,  fl]  =  compmodO 

clear 

clc 

k  =  3;  %  number  of  blades,  up  to  7 

syms  Phil  Phi 2  Phi 3  Phi 4  Phi 5  Phi 6  Phi 7  Omega  t  alpha 

tempPhi  =  [Phil  Phi 2  Phi 3  Phi 4  Phi 5  Phi 6  Phi7] ; 

Phi  =  tempPhi (1 :k) ; 
psi  =  Omega* t  +  Phi; 

rl  =  sym( ' rl ( t) ' ) ;  r2  =  sym(  ' r2 ( t )  ' ) ; 
r  =  [rl  r2 ]  ; 

Tlr  =  [100;  0  cos(rl)  sin(rl);  0  -sin(rl)  cos(rl)]; 

T2r  =  [cos(r2)  0  sin(r2);  010;  -sin(r2)  0  cos(r2)]; 
ml  =  (Tlr*T2r ) . ' 


for  i  =  l:k 

m2(:,:,i)  =  [cos(psi(i))  sin(psi(i))  0;  -sin(psi(i)) 
COS (psi (i) )  0 ;  0  0  1] ; 
end 

zetl  =  sym( ' zetl (t) ' ) ;  zet2  =  sym( ' zet2 ( t) ' ) ;  zet3  = 
sym ( ' zet3 ( t ) ' ) ;  zet4  =  sym( ' zet4 ( t ) ' ) ; 
zet5  =  sym( ' zet6 (t ) ' ) ;  zet6  =  sym( ' zet6 ( t) ' ) ;  zet7  = 
sym( ' zet7 (t) ' ) ; 

tempzet  =  [  zetl  zet2  zet3  zet4  zet5  zet6  zet7  ] ; 
zet  =  tempzet (1 :k) ; 

betl  =  sym( 'betl ( t) ' ) ;  bet2  =  sym( 'bet2 ( t ) ' ) ;  bet3  = 
sym( 'bet3 ( t ) ' ) ;  bet4  =  sym( 'bet4 (t) ' ) ; 
bet5  =  sym( 'bet 6 ( t ) ' ) ;  bet6  =  sym( 'bet6 ( t) ' ) ;  bet7  = 
sym( 'bet7 ( t)  '  )  ; 

tempbet  =  [  betl  bet2  bet3  bet4  bet5  bet6  bet7  ] ; 
bet  =  tempbet (1 :k) ; 


for  i  =  l:k 

T3z(:,:,i)  =  [cos(zet(i))  sin(zet(i))  0;  -sin(zet(i)) 
cos (zet (i) )  0;  0  0  1] ; 


75 


T2b(:,:,i)  =  [cos (bet (i) )  0  sin(bet(i));  010;- 
sin(bet(i))  0  cos (bet (i) ) ] ; 

T3p(:,:,i)  =  [cos(psi(i))  sin(psi(i))  0;  -sin(psi(i)) 
cos(psi(i))  0;  0  0  1] ; 
end 

for  i  =  l:k 

m3(:,:,i)  =  (T3z ( : , : , i) *T2b ( : , : , i ) ) . ' ; 

m4 ( : , : , i )  = 

(T3z (  :  ,  : , i ) *T2b ( : ,  : , i ) *T3p ( : ,  :  ,  i ) *Tlr*T2r )  . ' ; 
end 


syms  h  el  R  mbl  mb2  mb3  mb4  mb5  mb6  mb7 

tempmb  =  [mbl  mb2  mb3  mb4  mb5  mb6  mb7]; 
mb  =  tempmb ( 1 : k ) ; 
ul  =  sym(  'ul ( t) ' ) ; 
u2  =  sym ( ' u2 ( t ) ' ) ; 

rhoFI_I  =  [ul  u2  0] ; 
rhoHF  =[00  h] ; 

rhoBuH  =  [el  0  0] ;  %  changed  from  vector  form 

rhoPBd  =  [R  0  0]  ; 
rhoHF_I  =  (ml* rhoHF.  '  )  •  '  ; 

for  i  =  l:k 

rhoBuH_I ( : , : , i )  =  ml *m2 ( : , : , i ) *rhoBuH . ' ; 
rhoPBd_I ( : , : , i )  =  ml *m2 ( : , : , i ) *m3 ( : , : , i ) *rhoPBd . ' ; 
rho ( : , : , i )  = 

(rhoFI_I+rhoHF_I+rhoBuH_I ( : , : , i ) +rhoPBd_I ( : , : , i ) ) . ' ; 
V(:,:,i)  =  diff (rho ( : , : , i) , t) ; 

Vsqr (  :  , : , i )  =  V(l,:,i)A2  +  V(2,:,i)A2  +  V(3,:,i)~3; 
TBk(:,:,i)  =  1/2 *mb (i ) *Vsqr ( : , : , i ) ; 

end 


o„  9,  o,  9.  o,  o„  o„  9. 9,  a  0,  9.  o.  9,  9,  a  Q,  9,  O,  9,  9-  9.  9.  cu  a  9,  Q,  9,  o,  ao,Q,aQ,o,o,a  Ct  o,ao,ao,o,aQ,Q,o,o,o,o,o,aa  q,  q,  q,  q, 
75  to  o  o  o  o  o  o  To  o  o  75  c>  a  0  0  o  o  o  o  o  o  "o  0  0  "6  o  'o  o  'o  ”6  o  'o  "5  'o  'o  €  ”6  o  'o  'o  “6  'o  'o  “6  'o  'o  “6  "6  "o  'o  'o  ^  ^  'o  'o  'o  'o 


%  Potential  Energy  of  kth  blade  (UBk) 


9,9,9,9.9,9,9,9,9,9,9,9.9.9.9.9,9,^9,^9.0.9.9,90,9,9.0,00,0,0.0.0,9,0,90-9,0,0,0.0.9,0,0.0.0.0.5.0.0,0,0.0,0.0. 
h  h  h  n  t  o  h  h  ^  t  t  o  ^00  h  h  h  hooo  h  b  o  'o  h  "o  o  'o  'o  -6  ?  Id  '6  'o  ”6  "6  'o  'o  ^  "6  6  '6  '6  'o  "6  6  'o  '6  'o  'o  ^  1  'o 


o,  a  o,  a  o,  a  o„  o,  o.  o,  o,  o,  a  o,  o,  o„  o,  a  o,  a 
t  o  1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 


syms  Kel  Ke2  Ke3  Kdl  Kd2  Kd3  Ksl  Ks2  Ks3  z 

syms  kf 11  kfl2  kfl3  kfl4  kfl5  kfl6  kfl7  kf31  kf32  kf33  kf34 
kf 35  kf 36  kf 37 
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syms  kill  kll2  kll3  kll4  kll5  kll6  kll7  kl31  kl32  kl33  kl34 
kl35  kl36  kl37 

syms  Ke4  Ke5  Ke6  Ke7  Kd4  Kd5  Kd6  Kd7  Ks4  Ks5  Ks6  Ks7 
Ket  =  [Kel  Ke2  Ke3  Ke4  Ke5  Ke6  Ke7]; 

Ke  =  Ket (l:k) ; 

Kdt  =  [Kdl  Kd2  Kd3  Kd4  Kd5  Kd6  Kd7]; 

Kd  =  Kdt (l:k) ; 

Kst  =  [Ksl  Ks2  Ks3  Ks4  Ks5  Ks6  Ks7] ; 

Ks  =  Kst (l:k) ; 

kf It  =  [kfll  kf 12  kf 13  kfl4  kfl5  kfl6  kfl7] ; 
kfl  =  kf It (1 :k) ; 

kf 3 t  =  [kf31  kf 32  kf33  kf34  kf35  kf36  kf37] ; 
kf3  =  kf3t (1 :k) ; 

kilt  =  [kill  kll2  kll3  kll4  kll5  kll6  kll7]; 
kll  =  kilt (1:3) ; 

kl3t  =  [kl31  kl32  kl33  kl34  kl35  kl36  kl37] ; 
kl3  =  kl3t (1:3) ; 

for  i  =  1 : 3 

UBkl  =  1/2*  (bet  (i)/'2*kfl(i)+zet(i)/'2*kll(i)  )  ; 

UBk3  =  1/4* (bet (i ) ~4*kf 3 ( i ) +zet ( i) ^4*kl3 ( i ) ) ; 

UBk(i)  =  UBkl  +  UBk3 ;  %  +  UBk3 ;  not  required  for  simple 
model 
end 


£•  9c  %  9a  %  9a  9a  Sk  9-  £•  9-  9-  9c  9'  9-  9-  9*  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  9-  Q-  9-  Q-  Q.  a  9.  a  a  a  a  a  a  o.  q,  a  a,  a  q.  a  o,  g,  o,  c_  a  a  o„  a  a  a 

ootooootoooototf /o^  otttyo'o'o'6'6  o  o  o  o  c'o  o  o'6'6-6'o'6'o'6'6''6  o''6-6'6'6'6l)-6^''b^,6'6 


•6  "6  'o  "6  "do  -6  "6 '6  'o  75  '6  "oo  'o  -6  6 


%  Dissapative  Energy 


9.^^^o^9.o,9,^aaao.^aaaaao.ao.aato.aaao,ao,r,aatao,aQ,ataaaaao,aaaaaaaaao. 

'o'o'o'o'o'o'o'o^'o'6'6'o'6'6'6'6'6'6^'6'6'6^^'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'5^'6'5^b5ig^^b5^g^ 


9. 9, 9, 9. 9, 9. 9, 9. 9, 9,  o.  9,  a  a  a  a  q,  q,  9. 
15  o  c  "o  Id  o  c>  "o  o  o  "0  0  0  75  o  o  "o  "o  "o 


syms  DB1  DB2  DB3  Czetal  Czeta2  Czeta3  Vzetal  Vzeta2  Vzeta3 

syms  cfl  cf2  cf3  cf4  cf5  cf6  cf7 

syms  ell  cl2  cl3  cl4  cl5  cl6  cl7 

eft  =  [cfl  cf 2  cf 3  cf 4  cf 5  cf6  cf7]  ; 

cf  =  eft ( 1 :k) ; 

cl t  -  [ell  cl2  cl3  cl4  cl5  cl6  cl7]  ; 
cl  =  clt (1:3) ; 

Czeta  =  [Czetal  Czeta2  Czeta3] ;  %  convert  to  7  blades 
matrix 

Vzeta  =  [Vzetal  Vzeta2  Vzeta3]; 
for  i  =  1 : 3 

DBk(i)  =  l/2*diff (bet(i) , t)~2*cf (i)  + 
l/2*diff  (zet(i)  ,t)/"2*cl(i)  ; 
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end 


9999999999999999999999999999999999990.  9  o,  a  o.  9  9  9  q  q  o.  o_  q  o.  o.  9  o_  o_  o_  a  o,  q  o„ 
o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  '<}'6'6'o'b1>'o1)'6'6'6'6'o'ot>'o'o'6'6'o  o"6'8'&'b'6'b"o'6'o'o''o'6'6'6"o 


.999999999999999999 

j'o'b'o'o^t'o'o'O'o'o'o'o'o^'oo'o 


u  uA^/i  toojLUii  lui  r  uociayc 

%%%%%%%%%%%9c%9;%%9-qqqqq9q9qq9q999999q99999999999999999999999 
o  o  o  o  o  o  o  o  o  o  o'o'ott'b'0'ot^'o-o'o'o'6'6'6^  o^'6/o'6'5'6:o:o'6'o'6'6'6'6'6^'6'6'6:o'6?'6'6'6'6t^'6 

o,  9,  9  9.  9  9  9.  o.  aaaa  9. 9. 9. 9. 9.  q,  o„ 

00000000  O  O  O  O  '  O  6  O  O  'O  'O  'O 


9.  9  9.  9. 9, 9,  P,  9  9  9  9.  a  a  9  a  a  a  a  9.  o„  9  o_  o„  o„  o.  o_  o„  9  9  o„  9  9  9  o_  9. 9  c 
o  ot'o'o'o  o^'o'o  o'o'o'b'o'o'o  ot  o  o  o'bt>'b'6'6'b'6'6'6'b'6^i 

SSS'S'S'SSS'S'S'S'S-S'S'S'S'S-S'S' 

O^'O'O'O'O'O’O'O'O'O'O'O'O'O'O'O'O'O 

%  Kinetic  Energy  of  Fuselage  (TF) 

qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqc 

'o'o^^^^^xt^'0'0'6'0'6'6'0'0'0'6'6'6'6'6'6'6'6'6'6'6'6'6'6'0'6'6'6' 

9-  S'  S-  9'  S-  S'  S-  9-  9-  9"  S-  <?-  o,  9  o.  q  9.  q,  o, 
w  o  o  o  o  o  ^  o  o  'o  'o  "6  "o  "6  "o  "o 


syms 

K1 

K2 

c  = 

[Cl 

c2  ] 

v  = 

[vl 

v2] 

M  = 

[Ml 

M2] 

K  = 

[K1 

K2] 

112 


TFt  =  1/2* (diff (ul, t) A2*M(l)+diff (u2 , t)A2*M{2) ) ; 

TFr  =  1/2* (diff (rl, t) A2*I11  +  dif f (r2 , t) A2*I22- 
2  *dif f (rl, t) *diff (r2, t) *112) ; 

TF  =  TFt  +  TFr; 

S'S-asas^SS'O'S-S'S-S'S'S'S'S'SSS'SSS'SS'S'SS'S'S'SSSSSSSSS'SS'SSS'SSSSSSSS'S'SSSSSS 
o  o  o  o  o  o  o  o  o'o't)'ti'q'o'o'q'o'o'o'o'6'o'o'o'o'o'o'o'o;:o'6'6'o'6'o'6'6'6'6'6'6'6'6'6:6'6'6'6'6'6'6'6'6'6>5^'6^io 

9.  9. 9.  9. 9. 9. 9, 9, q o, 9,  90,0,9 9  9. S-  S' 
a  o  o  o  o  o  o  o  o  o  o  o  o  "q  o  o  ^  'o  'o 

%  Potential  Energy  of  Fuselage 

Sc^^ScScScScSc^S'S'S'S'S'S'S'SSS'S'S'S'S'SS'S'S'S'S'S'S'S'S'S-S'S'SqqqqaqQ-qaqqqqqaQ-qQ.qQ.qq 
o  o  o  o  o  o  o  o  o'o'o  o  x>o  o  o  x>o  t>ooo  06  o  6  o  o  o  'o  0  '6  o  o  0'S  o'6'6'6"6'6'6  o'6''o'o'6^'6'6'6'b'6"o'6"ou'6 

99999999999999999999 

'o'ot'o'ot^'o'o'o'6'6'0'6'6'6'6'6'6'6 

syms  KTl  KT2  KRl  KR2 


UFt  =  l/2*ulA2*KTl  + 
UFr  =  l/2*rlA2*KRl  + 
UF  =  UFt  +  UFr; 


l/2*u2' 

l/2*r2' 


2  *KT2 ; 
2*KR2; 


9-  S  9-  9  S'  S-  S-  S'  S-  9-  S-  S'  S'  S-  S'  S'  S'  S-  9  S-  9  S  q  9  q  9  q  0-  q  q  S-  o,  o,  0,  q,  o,  q,  q,  o,  o,  q,  q,  o,  o,  o,  o,  o,  q,  0,  o,  o.  o,  o,  o,  o,  0,  cl  cl 

-6  '0  '6  'o  o  060  ^600  ^  "6  0  0  ^6  0  'o  'o  ^  "6  'o  'o  '6  '6  "6  "6  “6 1  "6  '5  '6  '6  "6  ^  "6  "6  "6  '6  '6  '6  'o  "6  '6  1 1  '6  ?  1 ?  t  o  1 1  u  t5 1 

%%%%%%%%%%%%%%%%%%%% 

%  Dissapation  Energy  of  Fuselage  (DF) 

%S;S:§:S'S;S;S'99  9  9999  9- 999999  S' 99999999999  9  9999999999  9  9  999999  Q„qaaQ.q 

o  o  ^  o  o  ot>9>T>^5T5^5l5x>  o  ouuut^,o  6'o^  6'o'6'6'6'6'6'6'6'o'6'6'6'6'5'6'6'6'6'6''6'6/6^'6'6'6'6'6'6'6'6u'6 

9  a  a  9  9  a  a  9  9  a  a  9  a  9  a  a  a  a  o.  q 
'o'o'ol)'b'6l)'o'6'6'b-6'o'6'6'6  o^-S-o 

syms  CT1  CT2  CR1  CR2  VTl  VT2  VRl  VR2 


DFtv  =  l/2*diff (ul, t) A2*CT1  +  l/2*diff(u2,t) A2*CT2 ; 
DFrv  =  l/2*diff(rl,t) A2*CR1  +  1/2 *dif f (r2 , t ) A2 *CR2 ; 


78 


DFth  = 

l/2*diff (ul, t) ^2*abs (diff (ul, t) ) *VTl+l/2*dif f (u2, t) ^2*abs (d 
iff (u2, t) ) *VT2; 

DFrh  = 

l/2*diff (rl, t) A2*abs (diff (rl, t) ) *VRl+l/2*dif f (r2 , t ) A2*abs(d 
iff (r2 , t ) ) *VR2 ; 

DF  =  DFtv  +  DFrv  +  DFth  +  DFrh; 


o  o  o  o  o  o  o  o  o  t>  o  o  t>  o  ^3  ^>  o  t>  o  'o  '6  6  6 1>  t$  o  o  t>  o  -6  6  6  "6  '6 15 15 1>  'o  '6  'o  "6  "6  '6  "6  "6  'o  "6  "6  'o  '6  '6  %  %  %  %  %  % 

9,9,o,9.o.9.o,^o,9.Q,o.Q,o.aaQ.aao.Q. 

v'Q'uuvv^vuv'O'Ouu'O'Ou'Ov'O'O 

%  Aerodynamic  (Generalized  Aerodynamic  Forces) 

%%%%% 


^  v  u'oj.icj.  o<.  -i-  zj  viz  ncj,uu_yuauuL^  ruibcb  / 
aaQ,5.aaaaGaao,o,ao,ao,ag,ao,Q,Q,aaao,Q,Q.o,aaa.Q.o,o,o,Q,o,o,om.om  o  o  o  o  o  o  o  o  o  o 

Ol  0_  O.  CL  O-  O-  o.  cl  o,  o.  o.  cl  cl  cl  o_  cl 


%%%%%%%%%%%%%%%%%%%%% 


syms  thetl  thet2  thet3  thet4  thet5  thet6  thet7  cdO  mu 
lambda  rhol  a  c 


thett  =  [thetl  thet2  thet3  thet4  thet5  thet6  thet7] ; 
thet  =  thett (l:k);  %  do  for  7  blades 


Vair  =  [mu*Omega*R,  0,  Omega*lambda*R] ; 
for  i  =  l:k 

V_I_t ( : ,  : , i )  =  ( -V  (  :  ,  : , i ) +Vair 

V_Bd_t ( : , : , i)  =  simplify(m4 ( : , : , i) *V_I_t (:,;,i)).!; 


end 

for  i  =  l:k 
UR ( : , : , i ) 
UT (  :  ,  :  ,  i ) 
UP( : ,  :  ,i) 
UU(  :  ,  : ,i) 

end 

=  V_Bd_t ( 1 ,  : , i )  ; 

=  -V_Bd_t(2, :,i) ; 

=  - V_Bd_t ( 3 , : , i ) ; 

=  sqrt  (UP(  :  ,  :  ,  i) /V2+UT(  : 

for  i  =  l:k 

aoa ( : ,  : , i)  =  thet (i) -atan(UP( : ,  : , i) /UT( : ,  :  ,  i) )  ; 
dFbeta ( : , : , i )  = 

l/2*rhol*a*c* (aoa( : , : , i) *UU( : , : , i) *UT(  : , : , i) - 
cdO/a*UU ( : ,  : , i) *UP( : ,  :  ,  i)  )  ; 
dFzeta ( : , : , i )  =  - 

l/2*rhol*a*c* (aoa (i) *UU ( : , : , i) *UP( : , : , i) +cdO/a*UU ( : , : , i) *UT 
( :  ,  :  ,  i )  )  ; 

Mbeta_k( : , : , i)  =  0 . 7 *RA2 *dFbeta ( : , : , i ) *cos ( zet ( i ) ) ; 

Mzeta ( : , : , i )  =  0 . 7*R~2 *dFzeta ( : ,  : , i)  ; 

end 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

9,^9,9.9,9,9,9,aaaaao.o.aao,o.a 

'O'6'6'o'6'6'o'6<o'o'6'6'o^6'6'6'6'6'6^6 


o  o  o  o  o  o  o  o  o  o  o  o  o  o  o  t 
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%  Derivation  of  Equations  of  Motion  by  Lagrangian  Method 

9-  Sc  Sc  Sc  Sc  Sc  2-  9-  9-  9-  9"  9-  9^  S-  9-  9-  9*  9-  9-  9-  9-  9-  9-  9-  Q.  o„  o,  q,  a  a  a  q.  q,  a  q.  a  q.  a  q.  o.  o„  a  a  a  o.  a  a  a  cl.  ol  cl  q.  a  a  a  o  o  o 

o  o  o  o  6  -6  "6 1>  'o  *6  *6  *6  t>  -6  "o  "o  -6  'o  S  *6  -6  *b  *6  "o  -6  "o  "o  "o  "o  'o  'o  *6  "6  "6  *6  "6  *6  "6  *6  "6  "6  "6  'o  *6  %  %  %  yo  %  '6  '6  '6  %  %  %  %  %  %  % 


<9. 9, 9,  ct  o„  ct  o_  9,  a  o,  a  cl  9,  o.  ao.ae.ao, 
o  o  o  o  o  o  o  o'o'o  0  0  'o  'o'b'b'6'6'6 


syms 


DOFF  =  [ul,  u2,  rl,  r2]  ; 
DOFB  =  [  ]  ; 

DOFB  =  [bet,zet]; 


DOF  =  [DOFF,  DOFB]; 
dDOF  =  di f  f ( DOF ,  t )  ; 
ddDOF  =  diff (dDOF, t) ; 
[l,n]  =  size (DOF) ; 


ci,  o,  9. 9,  o,  a  o,  a  a  a.  a  9.  a  a  o,  o,  a  a  a  o,  q.  a  a  o,  a  o, 
o  o'b’D'ti  o'o^’o'o'o  o  o'6,6'6'6,'o'6'6'b'6'b 

9  9,  9, 9  o,  9.  9  9.  0.  q.  a  a  o.  o.  o.  o,  a  a  a  a 
o  o  o  o  o  o  o  o'o'o'oo  o'o'o 'o'o'o 'o'o'o 


'o'o'O'o'o'o'o'o'o'O'o' 


o  o  o  00  ooooooolc'b'o 


%  Transformation  between  time  dependent  and  independent 
notation  by  substituting  sets 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


,9.9,9,99.990.9999999990,99 

>ox>t>'o'0'o'0'0'oo6'6'o'o'6'6'o'o'6'6 


syms  ql  q2  q3  q4  q5  q6  q7  q8  q9  qlO  dql  dq2  dq3  dq4  dq5  dq6 
dq7  dq8  dq9  dqlO 

syms  ddql  ddq2  ddq3  ddq4  ddq5  ddq6  ddq7  ddq8  ddq9  ddqlO 
syms  qll  ql2  ql3  ql4  ql5  ql6  ql7  ql8  dqll  dql2  dql3  dql4 
dql 5  dql 6  dql7  dql8 

syms  ddqll  ddql 2  ddql 3  ddql 4  ddql 5  ddql 6  ddql 7  ddql 8 


DOFqt  =  [ql  q2  q3  q4  q5  q6  q7  q8  q9  qlO  qll  ql2  ql3  ql4  ql5 
ql6  ql7  ql8] ; 

DOFq  =  DOFqt (l:n); 

dDOFqt  =  [dql  dq2  dq3  dq4  dq5  dq6  dq7  dq8  dq9  dqlO  dqll 
dql2  dql 3  dql 4  dql 5  dql 6  dql 7  dql 8] ; 
dDOFq  =  dDOFqt (l:n); 

ddDOFqt  =  [ddql  ddq2  ddq3  ddq4  ddq5  ddq6  ddq7  ddq8  ddq9 
ddqlO  ddqll  ddql 2  ddql 3  ddql 4  ddql 5  ddql 6  ddql 7  ddql 8] ; 
ddDOFq  =  ddDOFqt (1 :n) ; 

setl  =  [  ]  ; 
set2  =  [ ] ; 
for  i=l:n 

setl  =  [setl 
setl  =  [setl 
setl  =  [setl 
set2  =  [set2 
set2  =  [set2 


%  size  of  DOF  vector 
maple ( ' ' = ' ' , DOF ( i ) , DOFq ( i ) )  ]; 
maple  (  '  '  =  w,  dDOF  (i)  ,  dDOFq (i)  )  ]; 
maple ( '  '  =  ' ' ,ddDOF(i) ,ddDOFq(i) )  ]  ; 

maple (''=’', DOFq ( i ) , DOF ( i ) )  ]; 

maple  ('  '  =  w  ,  dDOFq  (i)  ,  dDOF  (i)  )  ]; 
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set2  =  [set2  maple  (''  =  w,  ddDOFg  (i ),  ddDOF  ( i )  )  ]; 

end 

setl  =  maple (' convert setl set ') ; 
set2  =  maple (' convert ', set2, ' set ') ; 

T  =  TF; 

U  =  UF; 

D1  =  DF; 

GF  =  [0,0, 0,0] ; 

[l,n]  =  size (DOFB) ; 
for  i  =  l:k 

T  =  T+TBk ( : , : , i ) ; 

U  =  U+UBk ( i ) ; 

D1  =  Dl+DBk ( i ) ; 

GF  =  [GF  Mbeta_k ( : , : , i )  Mzeta (  :  ,  :  ,  i )  ]  ; 

end 

Temp  =  maple ( ' subs ' , setl , T) ;  %Temp  =  subs ( T , dDOF , dDOFq) ; 
save  for  time  comparison 

GFq  =  maple ( 'subs' , setl, GF) ;  %GFq  =  subs (GF , dDOF , dDOFq) ; 
save  for  time  comparison 

[l,n]  =  size (DOF ) ; 
for  i  =  l:n 

tempi  =  diff (Temp, dDOFq(i) ) ; 
temp 2  =  maple ( ' subs ' , set2 , tempi ) ; 
temp3  =  diff (temp2, t) ; 

LI  =  maple (' subs ', setl , temp3 ) ; 

L2  =  diff (Temp, DOFq(i) ) ; 
tempL3  =  maple (' subs ', setl , U) ; 

L3  =  diff (tempL3,DOFq(i) ) ;  %  check  dimensions  of  DOFq 
tempL4  =  maple (' subs ', setl , Dl) ; 

L4  =  diff (tempL4, dDOFq (i) ) ;  %  check  dimension  of  DOFq 
%  EOM(i)  =  simplify (Ll-L2+L3+L4-GFq(i) ) ; 

EOM(i)  =  (Ll-L2+L3+L4-GFq(i) ) ;  %  input  to  allow 

computation,  gags  on  simplify (EOM) 
end 


%  code  will  generate  EOM,  will  not  perform  operations  on 
EOM,  ie.  simplify,  eval,  etc. 

%  remainder  of  code  should  be  correct,  could  not  evaluate 
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9- 


o  Format  equations  of  motion  into  form  A  d2x/dt2  =  f 

2-S;9'9'9'9'9'$-$'$'9'^9-9'9-9'9'^9'9-9-9-9-9'9^^^aQ,o,o,o.o.o.o.o,o,Q,aQ,g.qag.nqao,Q.o,o,o,Q,o,o,Q,Q,o, 

'O'O'O'O'O'O^^'D'O'O'O'O'O'O'O'b'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O'O 


for  i  =  l:n 
for  j  =  l:n 

A (i ,  j )  =  maple ( ' coef f ' , EOM ( i ) , ddDOFq ( j ) ) ; 

end 


end 


[j,g]  =  size (ddDOFq) ; 
setZ  =  zeros(l,g); 

f(l:g)  =  -eval (subs (E0M(1 : g) , setZ , ddDOFq) ) ; 

Ax2dot  =  ddDOFq* A; 
for  i  =  l:n 

f(i)  =  (EOM (i ) -Ax2dot (i ) )  ; 
f ( i )  =  maple ( ' subs ' , sets , f ( i )  )  ; 
f  ( i )  =  - (simplify (f (i) ))  ; 

end 


9.9.9,9,9,9,q.ao,aaaaaQ,Q,aaaQ.aaaaaaaaaaa.aaQ.aQ,aaaaaaaaaaQ,aaQ.aQ.aQ,o,oJ,oJ,Q,oJ, 

^^^^^^'0'0'0'6^0'0'0'6^'0'0'0'0'0'0'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'6'0'^^^^'^'6'0'5'S 

% 


9. 9. 9,  a  a  a  a  a  a 
'O'q'6'6'0'6'6'6'6 


%  Change  notation  for  input  into  S- function 

cl  a  9-  9.  9-  9. 9-  9. 9,  o_  a  a  9.  a  a  a  9, 9,  9,  a  a  a  ao,aao,aaa  a  a  a  a  a  a  a  a  q,  cl  a  a  a  a  o„  a  a  a  o„  a  a  cl  o_  a  ol  cl  a  o.  cl 
"o  t>  "o  *6 t>  *6  -6  ^  -6  Id  t>  'b  'o  "o  t>  'o  t>  *6  *6  *6  *6  "0  t>  -6  'o  -6  *6  *6  "6  *6  *6  "6  ~o  t>  'S  'o  'o  "o  "o  'o  *6  *6  *6  "6  -6  -6  -6  "6  %  %  "6  "6  %  '6  %  -6  %  % 


syms  xl  x2  x3  x4  x5  x6  x7  x8  x9  xlO  xll  xl2  xl3  xl4  xl5  xl6 
xl7  xl8 


Xltemp  =  [xl  x2  x3  x4  x5  x6  x7  x8  x9  xlO  xll  xl2  xl3  xl4 
xl5  xl6  xl7  xl8] ; 

Xldot  =  Xltemp(l:n); 

Xl  =  Xltemp (n+1 : 2*n)  ; 

Xl  =  [Xltemp  Xldot] ; 

DOFqtemp  =  [DOFq  dDOFq] ; 

setX  =  [ ] ; 
for  i=l:2*n 

setX  =  [ setX  maple ( ' ' = ' ' , DOFqtemp ( i ) , Xl ( i ) )  ] ; 

end 

setX  =  maple ( ' convert ' , setX, ' set ' ) ; 

tempsetXl  =  [ abs ( 1 , xl )  abs ( 1 , x2 )  abs ( 1 , x3 )  abs ( 1 , x4 ) ] ; 


tempzeros  =  zeros ( size ( tempsetXl) ) ; 

A1  =  maple (' subs setX, A) ; 
fl  =  maple  (' subs setX,  f)  ,- 
f2  =  subs (fl, tempsetXl, tempzeros) ; 

syms  ul  u2  u3 
ut  =  [ul  u2  u3] ; 

Al  =  subs (Al,u,ut) ; 
f2  =  subs (f2,u,ut)  ; 
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o\°  o\°  o\° 


APPENDIX  F.  MATLAB®  FUNCTION  TO  ACCESS  THE  MAPLE  FUNCTION 
‘ABS’  TO  CALCULATE  THE  ABSOLUTE  VALUE 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
abs  .m 

This  function  is  designed  to  call  the  Maple  function 
'abs'  with  two  inpus t,  b  and  c  and  return  the  result 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  a  =  abs(b,c) 
a  =  maple ( 'abs' ,b, c) ; 
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APPENDIX  G.  MATLAB®  FUNCTION  TO  ACCESS  THE  MAPLE  FUNCTION 

‘SIGNUM’ 


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%  signum.m 

%  This  function  is  designed  to  call  the  Maple  function 
%  'signum'  evaluating  one  expression,  a 
%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 


function  s  =  signum (a) 
s  =  maple (' signum' , a) ; 


87 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


88 


o\°  o\°  o\°  o\° 


APPENDIX  H.  MATLAB®  FUNCTION  TO  ACCESS  THE  MAPLE  FUNCTION 
‘SIGNUM’,  TAKING  THE  DERIVATIVE 


o,  o,  a  a  9.  o,  a  9,  a  9,  g.  a  a  a  a  a  a  a  a  a  a  a  a  q.  o.  a  g.  a  a  a  g.  a  a  a  a  g,  a  a  a  a  a  a  a  o,  g,  a  a 
"o'6'6t>'6'bt>'c>'6-6'b-6'b'6'6'6'6'6'6-6'6'b'6'6'6'6'6'6'b*6'6'6"6'6-6'6'6'6'6'6'6'6'6'6'6'6-6-6'6'6'6'6 


signuml  .m 

This  function  is  design  to  call  the  Maple  function 
'signum'  using  a  variation  of  the  function  with 


a  =  1 

9.9,9.9,o,9.9,9.9.9,aaaaaq.aaaaaaaaaaao,aaaq.aq.aaaaaaaaaag.aaag.aaa 

^'6'o'o'6'6'o^-6-6'6‘o*6'o'o'6‘6'f5t>'0'6'o'6'6'6'6'6'6'6'6i5-6'6'6,‘6'6'6,'6'6'6'6'6't$'6'6''6,6'6'6'6'6'6 


function  s  =  signum(a,b) 
s  =  maple (' signum' , a, b)  ; 
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APPENDIX  I.  MATLAB®  FUNCTION  MUNION,  TAKING  THE  UNION  OF 
TWO  SETS  USING  THE  MAPLE  FUNCTION 

function  sym_union  =  munion (argl, arg2 ) 

stringa  =  char (argl);  %get  rid  of  extra  char's 
stringa  =  stringa ( 8 : end-2 ) ; 

if  length (stringa)  ==  0 
stringa  =  ' { } ' ; 

end 

stringb  =  char(arg2); 
stringb  =  stringb ( 8 : end-2 ) ; 

if  length (stringb)  ==  0 
stringb  =  ' { } ' ; 

end 


union_string  =  [stringa,'  union  ', stringb] ; 

union_string (union_string  ==  ' ] ' )  =  ' }  '  ; 
union_string (union_string  ==  '  [ ' )  =  '  {  '  ; 

union_string  =  [ 'maple ( ' '  ' ,  union_string,  '  '')']; 

tempi  =  sym(eval ( uni on_st ring) ) ; 

temp2  =  char (tempi) ;  % [tempi (tempi  ==  ' {  '  )  =  '  [ '  ] 

temp2(temp2  ==  '{')  = 
temp2(temp2  ==  '}')  = 

sym_union  =  sym(temp2) ; 
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